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1 . INTRQDUCTIQf^ AI)MD SUMMARY 

♦ 

Because of the possible launch configurations required to boost a space shuttle 
into orbit, it is anticipated that a large number of control effectors, including both 
aerodynamic surfaces and gimballed rocket engines, will be required to control the 
vehicle during ascent through the atmosphere. One objective in controlling the vehicle 
is to determine the deflection angle settings of the control effectors required to trim the 
vehicle for headwind and sidewind disturbances, and for bias torques due to solid rocket 
motor misalignments. Because of the launch configuration and the large number of controls, 
the control engineer is faced with two challenging problems. First, to compute the trim 
solution may entail solving a system of coupled, nonlinear equations. Second, if the 
number of control variables exceeds the number of independent trim equations to be 
satisfied, the trim solution is not unique. 

To solve the uniqueness problem, additional constraints must be imposed. A logical 
choice for the additional constraints is the minimization of a performance criterion that 
penalizes the degradation in vehicle performance caused by large trim deflection angles. 
The performance criterion used in this investigation penalizes the following effects: 

• Thrust loss (gain) by gimballing the engines away from their nominal condition. 

• Thrust loss due to drag caused by deflecting aerodynamic surfaces. 

e Excessive hinge moments on aerodynamic surfoces. 

o Large movement of the actuators for trim which hampers the flexibility 
needed for dynamic response. 

The inclusion of a performance criterion in the problem formulation results in an 
optimizotion problem with equality constraints to be solved for the trim solution. This 
formulation eliminates the uniqueness problem but the control engineer is still faced with" 
the problem of explicitly solving the equations for the trim solution. Furthermore, the 
control engineer is likely to want to perform the trim computations many times in order to 
consider changes in the following: 

© Flight regime (dynamic pressure) 

• Desired trim conditions 

® Launch vehicle configuration 
® Set of control effectors 

• Steady-state wind disturbances 

• Performance criterion. 


1 



To serve this need, a computer p^grarn entitled TRIMS was developed to solve the trim 
problem numerically. The equations for the trim solution are based on the method of Lagrange 
multipliers and in general are nonlinear. Two standard numerical methods, steepest-descent 
and Newton-Raphson, are available for solving the nonlinear equations. Application of 
these methods yields a pair of iterative algorithms for computing the trim solution that are 
included In the TRIMS program. The program user can select the desired method at the time 
of program execution. The Newton-Raphson method is more efficient for linear or nearly- 
I inear equations, but may fail to converge in severely nonlinear problems unless started 
near the optimum solution. If the trim equations are linear and the performance criterion 
is quadratic, then the trim problem can be solved explicitly. For this case the Newton- 
Raphson method converges to the exact solution in one Iteration. The current version of 
the TRIMS program for a Space Shuttle during ascent (described in Appendix C) solves 
the lateral trim problem. The lateral-directional dynamics are in the program and the 
required data (stability derivatives, moments of inertia, and etc.) supplied by MSFC are 
stored internally in the block data subroutine. The program permits multiple-case runs 
and the cost of computing a trim solution Is minimal. The program is in a modular form 
that facilitates changes in the data and/or the equations defining the trim problem. 

In computing the trim solution, the control engineer must specify the particular 
peiformance criterion to be used. There is no rule or theory for determining a unique 
performance criterion. The usual procedure is to vary the performance criterion and examine 
the different trim solutions that result. In the TRIMS program there are fourteen relative 
weighting factors in the input data that can be varied in the performance criterion. By 
varying these a family of acceptable trim solutions can be obtained for more detailed 
examination . 

Two methods for determining which of the acceptable trim solutions is preferable were 
considered. 

One possible method for selecting from among several trim solutions Is based on 
controllability. If the trim problem is nonlinear, then the controllability of the linear 
vehicle dynamics about trim will depend on the particular trim solution. In this case, the 
trim solution that results in the most controllable system could be used. The notion of a 
controllab ility Index Is developed. This index provides a criterion for ranking the trim 
solutions according to the degree of controllability. The controllability index is computed 
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from a symmetric, positive semi-definite contrp1|ldb^lity matrix and is defined as the ratio 
of the maximum eigenvalue to the minimum eigenvalue of the matrix. This ratio has a 
minimum value of unity for an orthogonal matrix. For an uncontrollable system, the 
controllability matrix is singular and the value of the controllability Index is infinite. 

A diffi cu Ity with using the controllability Index is that the controllability matrix is not 
unique and the value of the controllability index varies with the choice of this matrix. 

The controllability Grammain [I i is one possible choice for this matrix. Other controllability 
matrices are also considered in the development of the control lability index. A second 
method for selecting a trim solution is based on comparing the trim solution to the maximum 
allowable deflections. The rocket engines and aerodynamic control surfaces can only 
rotate a certain maximum angle. For particular flight times, a maximum hinge moment 
requirement can reduce the maximum deflection angle of an aerodynamic control surface 
below its physical limit. Obviously, the trim solution must be within these deflection 
limits. Moreover, to permit freedom of movement, a control deflection should not be 
too close to its angular limit. Hence, the requirement that all deflection angles be within 
their limits by a specified margin could be used to select the trim solution. For linear trim 
equations, a quadratic performance criterion with a diagonal weighting matrix can always be 
found for which the trim solution meets this requirement if such a solution exists at all. 

(This property of a diagonal weighting matrix might extend to nonlinear trim equations, 
but the more general case has not been studied.) The search for a trim solution satisfying 
this requirement can be accomplished by vorying the diagonal elements of the weighting 
matrix using the penalty function method discussed in Section 3.2. 


For the lateral trim problem of the Space Shuttle, there are two aerodynamic control 
surfaces (aileron and rudder) and five rocket engines (three orbiter engines and two solid 
rocket motors). The physical or hard limits on the aileron and rudder deflection are 


aileron 


I -40° 


rudder ± 30° 


and, as noted earlier, maximum allowable deflections can be less than these limits due to 
hinge moment restrictions which vary with flight time. The physical limits on the rocket 
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engine deflections were assumed to be 3 : 30>^ The maximum control deflections computed 
by the TRIMS program for the lateral trim problem exceeded the limits when the solid 
rocket motors are not gimballed. In numerous instances, a deflection angle exceeded 
100 degrees. When the trim constraint of zero net side force is removed, the maximum 
deflection angles decrease by an order of magnitude and are within the limits. 

In addition to the trim problem, the capability of the control system to damp out 
perturbations about trim must be considered. This can be identified as the dynamic response 
problem. In order to solve the dynamic response problem, it must be determined if the vehicle 
has sufficient dynamic control authority after trim conditions have been achieved. An 
approach to this problem based on the control labil ity Grammain used in defining the 
controllability index mention previously is studied. The controllability Grammain is used 
to compute the energy expended by each control in damping out errors from the trim 
conditions. This approach is limited, however, since it does not directly examine the peak 
deflection angles nor does it consider the realization of the feedback control system. 

It is advantageous to have a single method for solving both the trim problem and the 
dynamic response problem. The method should minimize the total control deflection 
tequired both to trim the vehicle and to damp out initial errors and random disturbances. 

If the vehicle dynamics are linear, optimum control theory provides the desired method. 

In Section3,4 the equations for the solution of the optimum control problem are derived 
for the case of bias inputs (trim problem) and random inputs (dynamic response problem). 

In Section 4.3, this theory is used to design an optimum feedback system for the lateral 
control of the Space Shuttle, and the closed-loop performance is simulated for a step 
change in side-slip angle. The computations for this example of the optimum control 
approach were performed with the aid of the Linear Systems Design (LSD) program 
developed at Singer-Kearfott under its Independent Research and Development program 
concurrent with this investigation. 

ft is recommended that further investigation of the trim problem for the Space Shuttle 
be performed with the aid of the TRIMS program for different combinations of controls (i.e. , 
gimbal solid rocket motors), performance criterion, and trim constraints. In addition, a 
more extensive design effort using the optimum control approach would merit further 
consideration. 
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2. PROBLEM DESCRIPTION 


The objective of this study Is to determine how the control effectors for the Space 
Shuttle con be optimally used to achieve trim and dynamic control In the presence of wind 
disturbances and bias torques due to misalignment of rocket engines. Launch vehicles have 
in the past been primarily controlled by gimballing the rocket engines. Various Space 
Shuttle configurations now under investigation Indicate that engine gimballing will not 
provide sufficient control to trim the vehicle for headwind and sidewind disturbances. 
Consequently, it may be necessary to use aerodynamic surfaces in conjunction with engine 
gimballing to achieve trim. ' Because of the severe cross-coupling problems encountered in 
the launch configurations. It appears that a large number of control effectors may be used. 

If the number of control effectors exceeds the number of quantities to be controlled, then 
the set of deflection angles to achieve trim is not unique. Thus, the control engineer in 
this case has a family (most likely on infinite set) of possible trim solutions to choose from . 
However, different trim solutions will result in different levels of performance and dynamic 
control. Consequently, the objective of the control engineer Is to select the trim solution 
that provides the highest level of performance and dynamic control . To achieve this a per- 
formance criterion, which ranks the trim solutions according to level of performance and 
dyrramic control, is defined. The problem then becomes, "What is the unique trim solution 
which optimizes the performance criterion?" 

The algebraic equations for computing the trim solution are derived from the differential 
equatioris describing the motion of the vehicle by substituting the desired trim conditions. If 
the number of control variables exceeds the number of (independent) algebraic equations, then 
the trim solution^is not unique. By addition of the performance criterion mentioned above, a 
meaningful optimization problem which can be solved for a unique trim solution, is obtained. 
This section develops the general problem In greater detail showing how the trim equations are 
derived from the equations of motion and the mathematical form of the perforrrrance criterion . 
The general equations for studying the dynamic response about trim are also derived. 

2.1 Trim Proble m 

In general the motion of the vehicle is governed by a set of nonlinear, time-varying, 
differential equations of the form 

X = a(x,t) + b( 6 , X , t) + c(x , z , t) + v(t) (2.1) 
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where x(f) = n x 1 vector defining state of the vehicle motion at time t 
6(t) = m X 1 vector of control deflections 
z(t) - X 1 vector of bias disturbances 
v(t) = n X 1 vector of random disturbances 
ci(x,t) = n X 1 vector function of x and t 
K 6/X,t) ~ n x 1 vector function of ^,x, and t 
c(x,z,t) ■ n X 1 vector disturbance function of x , z , and t. 

The trim problem is to find the set of control deflections 6^ that yield the desired steady 
state trim conditions x in the presence of bias disturbances . The bias disturbances 
model the effects of a steady wind and misalignment torques. The trim problem ignores the 
landom disturbances, i.e., v(t) 0 is assumed. Therefore, the trim solution must satisfy 

Let 

0 = a(x^,t) +E( 6^,x^,t) + c(x^, z^,t) (2.3) 

represent the subset of (2) required to calculate the trim deflections 6 where a" £ 

^ d ' ' 

and c are nxl vector functions with n<n . In orther words, in obtaining the algebraic 

equations in (2.2) from the differential equations in (2.1), it is possible that some of the 

equations in (2.2) are satisfied by x^ independent of 6^ . These equations, although 

used in (2. 1) for computing the dynamic response, are not used in computing 6^ and may 

be ehm mated from (2.2). This elimination which results in (2.3) replacing (2.2) will be 

illustrated by the lateral control of the Space Shuttle in Section 4. 

If m <n then no set of control deflections 6^ exists that satisfy (2.3), If m >n 
then a solution 6^ exists but is not unique. If m=n, there exist a unique solution , but 
unless b is a linear function of 6^ (i.e., b( 6 , x t) = B(x t) 6 . it may be difficult 
to find. 

For the case of infinitely mony possible trim solutions (m>£), certain solutions are 
preferable over others. An example of the latter Is a solution in which each deflection angle 
IS smoller in magnitude than for another trim solution. Trim solutions in which any of the 
deflection angles exceed the maximum allowable deflection should be excluded since such 
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solutions cannot be realized. Suppose additional constraints in the form of a performance 
criterion are included in the problem formulation. The solution that satisfies the trim 
conditions (2.3) and minimizes the performance criterion is unique. For this approach the 
trim design problem reduces to the appropriate selection of the performance criterion.. 

The performance criterion denoted by r is a scalar function of the control 

r = r(6j) (2.4) 

In the case m ^ n the trim problem is to find the set of control deflections 

= [ 6, , Sj , . . . , 6t^l 

that satisfy (2.3) and minimize the performance criterion (2.4), In general , (2 ,3) and 
(2.4) are nonlinear functions of 6^ and the resulting optimization problem with equality 
constraints can not be solved analytically. Numerical methods for solving the nonlinear 
trim problem are developed in Section 3,1. 

If the trim equation (2.3) is a linear function of ^ 

0 = ^(x^,t) +S(x^,t) +^(x^,z,t) (2.5) 

where B is a n by m matrix and if the performance criterion is a quadratic form 

r(6j) = 1/2(6^-6^)'R(6^-6^) (2.6) 

where r is a positive definite matrix and where 6 is the desired trim solution (in most 

' o 

instances 6 =0) then the trim problem is said to be linear , The linear trim problem can be 
o 

solved analytically and the equations are derived in Section 3.2. 

2.2 GENERAL CONTROL PROBLEM 

The trim problem is only part of the vehicle control problem. In addition to bias 
disturbances, the control system must be able to damp out sudden deviations from trim 
and to sustain proper vehicle motion in presence of fluctuating disturbances. An example 
is a sudden change or rapid fluctuation in the side wind velocity or equivalently the sideslip 
angle . The capability of the control system to handle rapid fluctuations in /8 , for 
example, is commonly determined by simulating the performance for a step change, impulsive 
change, or random noise with a specified frequency spectrum. The control system must be 
designed to maintain the control deflections within the physical limits and to return the vehicle 
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to trim within an acceptable setting time. This problem can be identified as the problem of 
dynamic response about trim. The first step in studying the dynamic response of the vehicle 
is to linearize the equations of motion about trim. Let Ax , A6 , Az denote deviations of 
the state, control deflections, and bias disturbances, respectively, from trim. 


Ax = X - x^ 
A6 = 6 - 

A z = z - z^ 


(2.7) 


Expanding in a Taylor series the nonlinear functions a, b, and c in (2. 1) about trim 
conditions results in the approximations 


a(x,t) - a(x .,t) + [ Ba/dx]Ax 
a 

b(x,t) = b(x^,0 + [ab/3x + [ ab/&6]A6 (2.8) 

Cl 

c(x,t) = c(x^,t) + C dc/Bx]Ax + [ dc/dzjAz 


Substracting (2.2) and (2.1) and substituting (2.7) and (2.8) yields the linearized equotions 
of motion 


Ax = AAx + BA 6 + C A z + v 


where 

A = Sa/Sx + db/dx + dc/dx (2.9) 

B = ab/36 (2.9) 

C = dc/^ 

Note that the partial derivatives are evaluated about the trim conditions and that for 

particular values of 6,#x,,z.,t the matrices A, B, and C are constant. 

d d d 

If the total motion (trim + dynamic response) is governed by I inear differential 
equations then (2.1) becomes 

X = Ax + B 6 + Cz + V (2.10) 

which has the same form as (2.9). The matrices A, B, C in (2.10) ore in general a function 
of time t . By considering only a number of fixed points along the trajectory the problem 
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reduces to a set of matrix equations of the form (2* 10) with constant coefficients. 

There are two general approaches for studying the general control problem including 
trim and dynamic response. 

Approach 1 1 First solve the trim problem for a set of acceptable trim solutions 
by varying the performance criterion (2.4). From this set select 
the particular trim solution that leads to the best dynamic response. 
Methods for determining the particular trim solution are developed 
in Section 3.3. 


Approach 2: Formulate a single performance criterion for the general control 
problem and solve for the optimum combination of trim solution 
plus dynamic response. This differs from the first approach in that 
two performance criteria are used in the former— one for the trim 
problem and one for the dynamic response problem. 

Consider all possible combinations of forces and moments that can be generated by 
the controls of the Space Shuttle. This set defines the control authority of the vehicle. 
The restrictions on the control authority are of the form of bounds on the deflection 
angle, i.e.. 


Tmin 


^ 6 .^ 6 . 

t tmax 


t = ^ f , . . f m (2.11) 


For most of the controls the maximum deflection and is the same in either direction 


I! i ^ 6 


tmax 



The primary problem is to find a control solution that satisfies the restrictions (2.11). The 
restrictions (2,11) are in terms of the total deflection angles resulting from both trim and 
dynamic response requirements. Hence, the second approach is preferable to the first 
approach. However, the second approach in general presents more difficult computation 
problems. If the equations governing the total vehicle motion are nonlinear then it may be 
necessary to use the first approach; the second may lead to an intractable problem. If, on 
the other hand, the equations for the total motion are linear, as is the case of the space 
shuttle dynamics in Section 4, then a design method in the category of the second approach 
results from the application of optimum control theory. The use of optimum control theory 
to solve the general control problem with both random and bias input disturbances is 
developed in Section 3.4 and the application to the lateral control of the Space Shuttle 
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is described in Section 4,3, 


Even when the control design is to be performed using optimum control theory, there 
are advantages to first solving the trim problem. The trim solution i$ much easier to compute, 
and sufficient control authority must exist to handle at least the trim problem. Furthermore, 
the solution to the trim problem can aid in the formulation of the optimum control problem. 
The correlation between the trim solution and the optimum control solution is considered 
in Section 3.4, 1 , 
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3* ANALYTICAL METHODS 


3. 1 ITERATIVE SOLUTION OF NONLINEAR TRIM PROBLEM 
3,K1 Lagrange Multipliers 

From (2*3) and (2.4) in Section 2 it was shown that the computation of the control 
deflections required to trim the vehicle for bias disturbances can be modeled as a problem 
of the following form: 

Find the vector 6 of dimension m which minimizes a scalar function of 6 

min r(6) 

6 

subject to a set of n equality constraints 

0 = a + b(6) (3.2) 

For simplicity, the subscript "d" has been dropped from 6^ and (2 .3) has been rewritten 
as (3.2) where 

a s a (x^ , t) + c(x^ , t) 
b = £'( 6, , t) 

n = n 

For a particular point in time t along the trajectory and for a particular set of desired 

trim conditions x and bias disturbances z , the vector a in (3.2) is a constant and 
d 

the vector b ia a function of 6 only. 

In order to achieve a well-defined optimization problem the performance criterion 
r(6) is assumed to have the following properties; assume that r is differentiable and 
let 6* be the value of 6 that minimizes r( 6) . (Here the subscript d has been 
dropped from 6^ since just the properties of the performance criterion r are of interest 
irrespective of the trim equation (3.2).) 

r(6)^r(6*)^0 (3.3) 


(3J) 


n 



Then in some neighborhoods of 6* / tbe performance ciferion has the property that the 
grociient srHsfies 

r = 0 for a = 6* 

Sr / d 6 S (3.4) 

[ 5 ^ 0 for 6 6* 

where ^r / S6 = [ Sr / S6, / - • • / 1 • Furthermore the second partial 

1 m 

derivative of the performance criterion or Hessian matrix satisfies ■} 

. « [ > 0 for 6 = 6* 

S r/S6^ ( . (3.5) 

[ 0 for 6 6* 

where ( S^r / S6^). . - / S6 j S6 . 

The basic approach for solving the nonlinear trim problem given by (3.1) and (3.3 is to 
apply the well-known method of Lagrange multipliers. 

Define a new scalar function h (the Hamiltonian) by 

h(6, X) = r(a) + V(a+b(6)) (3.6) 

where X is a vector of n unknown parameters, commonly referred to as the "Lagrange 
multipliers". The fundamental idea underlying the method of Lagrange multipliers is that 
if 6* / X* is the solution that minimizes h then 5* is the solution that minimizes r 
and satisfies (3.2). 

Assuming that the functions r(6) and b(6) are differentiable, the equations for 
the minimal solution 6* , X* can be obtained by differentiating h and setting the 
derivatives to zero. This gives 

dr / + X^ 6b/ S 6 = 0 (3.7) 

a + b( 6) =0 (3.8) 

-|- The notation " >0" means the matrix is positive definite and " ^0" means the matrix is 
positive semi-definite. For reference purposes see Appendix A for a discussion of 
differentiation by a vector. 
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This is a system of m+n equations in m4n unknown 6 and X . Only in special cases 
can (3.7) and (3.8) be solved explicitly. In general, numerical methods must be used to 
solve (3.7) and (3.8), Iterative numerical methods for determining the solution 6* , 

X* that minimizes h given by (3.6), start with an initial guess 6^ , X^ and then 
proceed to compute a sequence of solutions 


X^ t ^2 / • • »/ X|^ / • • • 


which converge to the exact solution 6* , X* 




Two such numerical methods are described in Sections 3.1.2 and 3.1.3. 

3.1.2 Numerical Solution by Steepest Descent Method 

One numerical method. In common use for many years, for finding the minimum of a 
function is that of "steepest descent". The steepest descent method is a 1st order gradient 
method and uses an iterative algorithm for improving the estimate of the solution so as to 
come closer to satisfying the zero slope conditions 

dh/36 = 0 and dh/ SX = 0 

The method computes value of Xj^ is not used to continue 

the iteration. The method partitions the vector 6 according to 


a' = [x lu ] 

^ 

n m-n 


where the subvectors x and u are computed separately. 

Application of the steepest descent method gives the following steps for computing 
\+i ' “k+i \ ' “k • 
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1) From X|^ , Uj^ compute the column vector b( 6) . 

2) From Xj^ , Uj^ compute the matrices Sb / 3x , / Su . 

3) Compute the new estimate of subvector x according to 

Ax, =-(3b/3x) \a +b(6)) 


\+i = \ 

4) From t compute the row vectors 3r/3x , 3r / 3u * 

5) Compute the vector of Lagrange muitipUers according to 

^k+l ~ ^b/3x) ^ 

6) Compute the gradient of h with respect to u using 

3h/du =3r/3u + Xj^^^(3b/3u) 

7) Compute the new estimate of subvector u according to 

A Uj^ = -{j(3h/3u)' 

8) Repeat steps 1) through 8) with the updated solution Xj^^^ , u^^^^ until 
the total error is very small 

lUxu 11^ + 11 n ^ 

where the norms are given by 


II II 

i AU|^ 11^ = 


A flow chart showing the basic steps required to implement the steepest descent 
method for solving the trim control problem on the computer is given in Figure 3.1. A 
graphical Interpretation of first order gradient methods is given on p. 20 of [ 2 ]. 
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Firsf order gradient methods usually show substantial improvements in the first few 
iterations but have poor convergence characteristics as the optimal solution is approached. 
A second— order gradient method/ which uses the "curvature" as well as the . slope at. the-, 
nominal point/ is discussed in the next section. Second order gradient methods have 
excellent convergence characteristics as the optimal solution is approached but unless the 
initial guess is in the region of convergence then the method may not converge or may 
converge to the wrong solution. 

3.1.3 Numerical Solution by Newton-Raphson Method 

Newton— Raphson method (or second— order gradient method) for locating the minimum 
point of a function uses both the first and second derivative at the nomirxal point to extrap- 
olate a new estimate of the solution. A detailed description of the Newton-Raphson 
method is given in [ 1 ]. 


Using the Newton-Raphson method to find the minimum solution of h( 6 / X) given 
by (3.6) yields an iterative algorithm for computing the trim solution. To obtain the 
equations for computing from f / first expand h( 6 , X) in a Taylor 

series about , X|^ . . 


h(6 / X) = h(6 / X|^) + 







^66' '’6X 


h.i h 


A6 

.,1 

A6 

f 


A6 




^ m 




— — - 

0, X 


AX 

AX 


^xsl '’xx 


AX 


+ . 


(3.9) 


where 


A6 = 6 - 6, 


(3.10) 


AX = X -X 


Differentiating (3.6) gives the following set of equations for evaluating the derivatives 
In (3,9) 


h^^ dh/36= Sr/66-tX'Ob/36) 

0 

h = 3h/6X =a' -tb'(6) ( 3 . 11 ) 
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hjg = = 3^r/S6^ +X'(a^b/36^) ' 

h =3^h/&6aX = (6b/a6)' 

h. J = 6^h/axa6= ab/d6 (3.12) 

h = a?h/ax^ = 0 

From (3.9) the equations for computing the new estimate of the solution are: 

=6|<+Aa|< 

(3.13) 

^k+l “ 

where the incremental corrections A6|^ / AXj^ are the solution of a system of linear 
equations 


‘’sa ; '"ex 


1 

, > 

' 

1 

J 

h. . , h 


A A, 

xa ' XX 

L. 1 J 


k 



(3.14) 


Note that the derivatives are evaluated about the nominal point • 


To summarize, the steps in the Newton-Raphson method for computing 


'k+r 


from 


®k ' 


are as follows: 


1) From 6j^ compute the column vector b( 6) . 

2) From 6, compute the matrix Sb/B6 . 

^ 2 2 

3) From 6|^ compute the tensor B b/d6 . 

4) From 6j^ compute the row vector Br/d6 . 

2 2 

5) From 6^^ compute the symmetric matrix B r/B6 

6) Compute the first-order gradient terms h , h according to (3. 11). 

V A. 

7) Compute the second-order gradient terms h „ , h , h according to 

05 M AO 

(3.12). (Note that h = 0.) 

AV 

8) Compute the incremental correction to the solution by solving (3. 14) which gives 


) 
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(3.J5) 


A6. = - ER‘’-R‘’B'(BR"’Br’BR‘’]h^- [R"’B'(BR'Vr’ 

= - [(BR"’b')"’br"’ ] h'g+[(BR‘’B')‘’]h' 


where the matrices R and B are defined by 
" ^86 

9) Update the solution according to (3.13). 

10) Estimate the error In the solution by computing the norms 


I f = * 6k 

II 


(3.16) 


11) Repeat steps 1) through 11) with the updated solution until 

the sum of the norms is very small as given by 

II 11^ + iAX,^ 11^ <€ 

A flowchart showing the basic steps required to implement the Newton-Raphson 
method on the computer is given in Figure 3.2. 
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Fiaure 3-2 Flow Chart of Newton -Raohson Method for Solving the Trim Control Problem 
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3.2 SOLUTION OF LINEAR TRIM PROBLEM 
3.2.1 Explicit Formulas 

In the previous section the general nonlinear trim problem defined by (3. 1) and (3.2) 
was discussed. The case of a linear trim equation 

a+B6 = 0 (3.17) 

and a quadratic performance criterion 

r(6)= 1/2(6- 6^)'R(6-a^) . (3.18) 

is referred to as the linear trim problem and can be solved explicitly. The scalar Hamiltonian 
function corresponding to (3. 17) and (3.18) Is 

h(6A) = l/2(6-a )'R(6-6 ) +V(a+B6) (3.19) 

o o 


The vectors and matrices in the right-hand side of (3.19) are defined below 

5 = m - vector of control deflections. 

5 = m - vector of desired control deflections, 
o 

X = vector of Lagrange multipliers of dimension (m-n). 
a = constant vector of dimension n . 

B = constant matrix of dimension nxm . 

R = constant positive definite matrix of dimension mxm , 


The trim solution Is computed by determining the values of 6 and X that minimize 
the scalar function h . Differentiating (3.19) and setting the derivatives to zero gives 


(Bh/d6)' = R(6-6 ) +B'x = 0 
o 

(bh/dX)' = a + B6 


(3.20) 


The vector-matrix form of (3.20) is 



(3.21) 


Premultiplying both sides of (3.21) by the inverse of the square matrix on the right hand side 
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(3.22) 


of (3,21) gives that the optimum trim solution is 

6 = ]6 - B^a 

where 

B# = R"’b'(BR”’b')"’ (3.23) 

Note that B^ is a right inverse of B (i.e., BB^ = t) . Substituting (3,22) and (3.23) into 
(3. 18) and (3. 19) gives that the minimum values of performance criterion and Hamiltonian 
function are 

h = r = l/2(o+B6 )'(BR‘’B')'\a+Bfi ) 
o o 

Consider the example of triming sidewind induced roll and yaw moments using aileron, 

rudder, and the yaw deflection of a single rocket engine. Setting the rolling and yawing 

moment coefficients to zero (C. = C =0) gives in vector form 

n 



The trim equations given by (3.24) or (3.25) are a set of 2 linear equations in three 
unknowns there is one more unknown than equations, (3.24) has 

an infinite family of possible trim solutions. 

A graphical representotion of the possible trim solutions can be seen by depicting 
(3.24) in the yaw-roll moment coefficient plane as shown in Figure 3.3. The four vectors 
formed by the stability derivatives are represented by solid arrows where the following 
numerical values were chosen for the example 
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Sy 


- 0.20 . 

0.80 


c 


0.10 

IR 



C 


-0.10 

nR 




1 

n 

1 

1 

0.07 

1 

o 

1 


0.08 


1 — 

u 

u 


-0.67 



0.35 

, * 


The dotted cupve in Figure 3.3 represents one of the trim solutions for the case, ^ = 1® , 
and is the vector diagram corresponding to the left hand side of (3.24). The values of the 
deflection angles are 


6^Y = -^^'5 = 2.807 6^ = 4.133 

and are equal to the lengths of the dotted arrows divided by the lengths of the corresponding 
(parallel) solid arrows. The sign of the deflection angle is positive if the dotted arrow and 
the corresponding solid arrow point in the same direction, and the sign is negative if the 
directions are opposite. 


The addition of a performance criterion to be minimized will yield a unique solution 
for (3.24). For illustration, one possible choice might be 


r( 6) = ( +( 6 „/ 20°)2 +( 8 ^/, 0«)2 


(3.26) 


where 15® , 20® , 

. and 10® 

are the corresponding maximum deflections. 

From (3.25) it 

follows that for j3 

= 1® 





B = 

-0.20 

0.10 

0.07 


(3.27) 


0.80 

-0.10 

0.08 

0.35 



L 


- 



and from (3.26) 



- 




1/225 

0 

0 



R = 

0 

1/400 

0 


(3.28) 


0 

0 

1/100 




Substituting (3.27) and (3,28) into (3.22) and (3.23) gives the solution 



®EY 


0.10® 

6 = 

«R 


5.70® 


«A _ 


1.72® 


(3.29) 
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As an example of how changes in the performance criterion effect the minimal solution 
suppose in place of (3.26) 

r(6) = (6^y/20°)^ +(6p/20°) +(6^/20°)^ 

then 
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3,2,2 Performance Criterion Selection 

When infinitely mciny trim solutions are possible, certoin solutions definitely require 
more control authority than other solutions and should not be used. In particular, given 
a trim solution 6 , if it is possible to find another trim solution 6 * such that for each 
control 

1 I < I 6, 1 t = 1 , 2 , . . . , m (3.31) 

where the strict inequality holds for some controls then 6 should not be used. Property 
(3.31) partitions the possible trim solutions into. two disjoint sets. If 6 satisfies (3.31) it 
will be referred to as an unfavorable trim solution and if 6 does not satisfy (3.31) it will 
be referred to as a favorable trim solution. The problem of selecting a form of the performance 
criterion that guarantees a favorable trim solution hos been solved. 

At this point a simple example is helpful in studying the properties of the trim 
problem. Suppose there is a single trim equation 

0 = 6-26^+62 ( 3 ‘ 32 ) 

with two controls 6 ^ and 62 • The general form of the performance criterion for the case 
of two controls is 

r = l/^r,6^ + l/2r262 +r3 6, 62 ( 3 . 33 ) 

where 



and 

r^ > 0 , •'2 ^ ^ 

One approach for graphically representing the trim problem is to consider 6 - [ 6 ^ , . . . , 6 ^ 3 ^ 
as defining the coordinates of a point in an m-dimensional space which sholl be referred to 
as the solution spoce . This approach is different from the graphical representation in 
Figure 3,3 where each coordinote corresponds to one of the scalar trim equations and hence 
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might be referred to as the equation space representation. For this example the loci of 
possible trim solutions in the solution space is the straight line defined by (3.32) and shown 
in Figure 3.4. The segment of the straight line between points P and Q defines the set 
of favorable trim solutions and the remaining two segments on either side of P and Q 
define the set of unfavorable trim solutions. 

For each fixed value of the performance criterion, there corresponds a closed contour 
cut^e in the solution space. For (3.33), r = constant defines an ellipse centered at the 
origin of the solution space. By parametrically increasing the value of r a family of 
concentric ellipses of increasing size Is generated. One of these ellipses will be tangent 
to the straight line passing through P and Q , The point of tangency is the optimum 
solution. For example suppose " '^2 ~ ^ *^3 ^ ^ then the loci of constant performance 

are circles as Illustrated in Figure 3.4(a) for r - 0.5 and r - 3.6 . The circle with 
r ^ 3.6 intersects the straight line at the single point 6^ =2.4 and ^ - 1 .2 . 

This is also the optimum solution obtained using the formulas (3.22) and (3.23). For the 
case r^ = 4, '*2 ^ ^ ^ ^3 ^ ^ the optimum el I ipse Is 

2 2 
18 ^46^ + 62 

and is tangent to the straight line PQ at 6^ = 1.5 and ~ *3*0 . 

The above example illustrates how varying the weighting matrix R in the performance 
criterion leads to different trim solutions. However, there are more ways of varying R 
(degrees of freedom) than necessary. This means different choices of the R matrix can 
lead to the same optimum trim solution. 

The redundancy in the selection of R suggests that R can be restricted to a diagonal 

matrix without disregarding a favorable trim solution. This assumption simplifies the selection 

of R For the example illustrated in Figure 3.4, the principle axes of the ellipse will 

coincide with the coordinate axes in the solution space when and only when R is diagonal 

(i.e,, = 0) . As a consequence, the optimum trim solution for a diagonal R matrix will 

o 

always lie on the line segment PQ (region of favorable trim solutions). The optimum solution 
point in Figure 3.4 will move from point P to point Q as the ratio of the diagonal elements 
r^ / r 2 increases from 0 to ■«'’ . Thus Increasing the weighting on 6-| relative to the 
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weight-ing on 62 causes 1 6 ^ 1 to decrease and | 62 | to increase. 

This example illustrates the following general properties of the weighting matrix 
in the performance criterion; . 

Property 1; The optimum trim solution fora diagonal R matrix is always 
g favorable trim solution. 

Property 2; Any favorable trim solution is the optimum solution for some diagonal 
R matrix. 

The general proof of the first property is hot difficult. Let 6 be the optimum 
trim solution for 

R = Diag [r^ , . . . , r^] 

then the minimum value of the performance criterion is 

r = 6, + . . . + r 6^ ) (3.34) 

II mm 

Suppose 6 is an unfavorable trim solution, then there exists another trim solution 6 ^ 
satisfying (3.31) for which the value of the performance criterion is 

r* = + . . . +r 6*f ) (3.35) 

II mm 

Comparing (3.34) to (3.35) term by term, it follows from (3.31) and > 0 that 
r* < r 

But r is the minimum value and hence a contradiction! Therefore, 6 cannot be an 
unfavorable trim solution. 

Given a favorable trim solution 6 it should be possible to find a diagonal R matrix 
in the performance criterion for which the optimum solution is 6 . A general method for 
constructing such an R matrix or equivolently, a general proof of the second property has 
not yet been found. 

The formulation of trim control problem given by(3. 17^ and (3.18) is an optimization 
problem with equality constraints. However, as pointed out in Section 2, inequality constraints 
also exist due to the physical limitation on the control deflections. For a symmetric control, 
these will have the form 
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I 6. I ^ 6. t = 1, . . . , m (3.36) 

' t ' imax 

6 = maximum allowoble deflecHon of the ffh control 

t max 

The Inequality constraints are not included explicitly in the problem formulation since an 
optimization problem with both equality and inequality constraints is difficult to solve. 

Instead, the inequality constraints are handled by the penalty function method.. 

The basic idea of the penalty function method is to repeat the computation of the 
optimum trim folution for different R matrices in the performance criterion until each ratio 
16.1/6 'S l©ss than one and the difference 6. - | 6, 1 is sufficiently large to 

provide the additional control required to solve the dynamic response problem. The 
procedure for varying the elements of R is simplified if R is restricted to be a diagonal 
matrix. From the properties of a diagonal R matrix discussed previously, this restriction 
does not exclude any favorable trim solutions but does exclude all unfavorable trim solutions. 
As an illustration of how to vary the diagonal elements of R , suppose the pptimurri . 
solution for ' 

R = DiagCr,, ♦ ♦ ♦ / r 1 
I m 

results in one of the deflections 6. exceeding its limits. The next step is to increase the 

t/ 

corresponding weighting factor r. and solve the problem again. Repeat this procedure 

until 6 . is smaller than the maximum deflection. An increase in the weighting factor r 
t t 

will cause the magnitude of 6 to decreose at the expense of increosing the magnitude 

C/ 

of other deflection angles. If no adjustment of the weighting factore results in all the control 
deflections being within their corresponding limits then the launch configuration does not 
possess sufficient control authority. If the limits are exceeded for every control then from 
Properties 1 and 2, mentioned earlier, no acceptable trim solution exists. 

The modification of the performance criterion to produce a more desirable trim solution 
can be focilitoted by realizing that for small perturbations the change in the optimum trim 
solution is proportional to the change in the weighting factors of the performance criterion. 
Computing the differentials of (3.22) and (3,23) for the case 6^ = 0 gives 

d6=-dB^.a (3.37) 

-dB# = (j-B#B)R"’- dR • B# (3.38) 
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The derivation of (3.38) makes use of the identity 
d(R"’) = - R"' • dR • R"’ 

From (3.23) and (3.38) it can be shown that 
B ' d6 = 0 

which also follows from computing the differential of (3.17), Equations (3.37) and (3.38) 
showthat for small perturbations d6 varies linearly with 6R . Let d6^ denote the change 
in the trim solution due to dR^ , i.e., 

dR^-d6, 

Substituting 

dR = 

into (a 24) where w. Is an arbitrary scalar results in 
dfi = 

Thus, replacing R by R + dR causes the optimum trim solution to become 6 + d6 . 
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3.3 CONTROLLABILITY AND DYNAMIC RESPONSE 
3.3.1 Contro I lab M i fy Gramm tan 

The use of the controllabUity Grammian for studying dynamic response about trim 
is developed below. The trim solution uses part of the control authority. If the vehicle 
deviates from trim due to random disturbance or a sudden wind gust then it must be 
determined if the control effectors have sufficient authority in reserve to return the vehicle 
to trim. By using a different trim solution, better dynamic response performance could 
possibly be ochieved with respect to the control limits. The problem of determining which 
controls are most effective in zeroing out deviations is also of interest. If there are more 
control effectors ovailable than required it may be possible to disregard those controls 
whose effectiveness is small. 

Basic Theory 

In vector-matrix notation the linearized equations of motion about trim hove the general 

form 

X ~ Ax + Bu (3.39) 

where 


X = state vector of dirnension n 
u = control vector of dimension m 

The equation for the solution is 

t 

x(t) = <^(t) x(0) + r 4>(t-T)Bu(T)dr (3.40) 

b 

where the transition matrix is 

$(t) = (3.41) 

The control signal that will drive the error to zero at time T is 

u(t) = . B'*'(-t)W"’x(0) (3.43) 

where 

T 

W s W(T) = J ^>(-t)BB'^'(-t)dt (3.43) 

0 
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The matrix function W(t) Is referred to os the "controllability Grammlan" [11 
Substituting (3.42) Into (3.40) and using (3.43), It can be shown that x(T) = 0 

A useful criterion for indicating the amount of control effort Is given by the integral 


E = r u'udt 

0 


(3.44) 


which may be viewed as proportional to the total "energy" expended by the control effectors 
in returning the vehicle to trim. Substituting (3.42) into (3.44) and using (3,43) yields the result 




E = x'(0)W x(0) 

Thus the controllability Grammlan W(t) provides a means for computing E , 
Let denote the "energy" expended by the tth control effector, then 


(3.45) 


E, = f u?dt 
i Jq i 


t = 1. , . m 


(3.46) 


where 


u^(0 = B'$'(-0W"’x(0) 


and is the tth column of the B matrix. Substituting (3.47) into (3.46) results in 


with 


= x'(0)W"'w,W''x(0) 


Wj = J" 


(3.47) 


(3.48) 


(3.49) 


The ratio E^/E is a convenient measure for determining the relative effectiveness of the 
tth control effector. Upon substituting 


BB' = B, B'. + B„B' + . . . + B B' 3.50) 

112 2 mm 

into (3.43), it follows from (3.45), (3,48), and (3.49) that 

W = W, +W., + , , . + W (3.51) 

1 2 m 
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and 


E - Ej + + . . . + 


(3.52) 


Another opproach for computing E and E^ is obtained by rewriting (3.44) as 


E = trace < f ou'dt 

0 

Substituting (3.42) into (3.53) gives 

E = trace { B'MB} 


where 


(3.53) 


(3.54) 


-1 -1 

M= J $'{-t)W x(0)x'(0)W '^(-t)dt 
0 

Repeating this approach for (3.46) and (3,47) leads to 


(3.55) 


. E^ = = 1, 2, . . . , m (3.56) 

The advantage of using (3.56) in place of (3.49) is that instead of computing W,, Wj, . • ; , 
only have to compute M . The disadvantage is that if the initial state vector x(0) changes 
then M must be recomputed where os the matrices are not a function of x(0) and hence 
do not change. 

Computation of Controllability Grammian 

Several methods for computing the matrix W = W(T) defined by (3.43) are discussed below, 
Eigenvector Transformation 

Suppose a new set of state variables q(t) are introduced that ore related to x(t) by 

q=Qx (3.57) 

where by assumption Q is a nonsingular matrix. 

Substituting (19) into (1) gives 

q = + Bu (3*58) 
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where 


A = QAQ 


B = QB 

Let W(t) denote the controHabMity Grammian computed from (3.58) then defining 
^ = '^(T) ond applying the definition (3.43) to (3.58) results in 

W = QWQ' or W = q"’wq“ (3.59) 


If 


A = Diag L\y X^, - . . / ] 

where X^ are the eigenvalues of A then the columns of Q ^ form the corresponding set 
of eigenvectors. In this development it is assumed that the eigenvolues are real and distinct. 
The method can still be applied to the complex and the multiple eigenvalue case but the 
computations are more complicated. This method will not be generalized because it is 
intended only for illustration purposes and as a means for checking the other methods. If 
A is a diagonal matrix then the tronsition matrix is a diagonal matrix with diagonal 
elements. 




t = 1,. . . , n 


which upon substitution into the definition of the controllability Grammian (3.43) gives that 
the element of matrix W in row t and column J is 


w, ,= b'b, J- 

where b. is the vector of dimension m formed by the tth row of B, i.e., 

= 


(3,60) 
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Integrating (3.60) gives 




-(X^ + X_^)T.,]/_^^ ) _ ( 3 . 61 ) 

J ^ * • *f ^ 


Combining (3.59) and (3.61) defines the eigenvector transformation method for computing W 
To illustrate, consider the exomple 



-1 

2 

0 


1.5 

-1 

A = 

0 

-3 

0 

B = 

-4 

2 


1 

2 

0 


-10.5 

-25 


T = 0.5 
m = 2 
n = 3 


If the transformation matrix and its inverse are 


Q = 


then 


-2 -2 0 



-0.5 .0 

0.25 0 -0.25 

Q = 

0 0 

0 1 0 



-0.5 -4 

L 

J 



-1 0 0 



5 -2 

0 -2 0 


B = 

3 6 - 

0 0 -3 



-4 2 

-- 


A = 


Substituting into (3.61) 


W = 


W = 


Numerical Integration 


24.92 

3.48 

-38.33 

3.48 

71.88 

0.0 

-38.33 

0.0 

63.62 

W into (3, 

59) 


31.51 

-44.45 

38.48 

-44.45 

63.62 

-44.45 

38.48 

-44.45 

1 195.47 


-1 


-1 


(3.62) 


(3.63) 


Let F(t) represent the integrand of (3.43), I.e., 
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(3,64) 


T 

W = J Fdt 
0 

Differentiating F results in the matrix differential equation 
. F = AF + FA' 

Integrating both sides of the above equation from 0 to t gives the following linear 
matrix differential equation for computing the controllability Grammian W(t) 

- W = AW + WA' - BB' , W(0) =0 ' ' (3.65) 

The solution to (3.65) at t = T is the value of the integral (3.43). Similarly, the linear 
matrix differential equation for computing W^(t) is 

-W^ = AW^ +W^A' - B^B'^ , W^(0) =0 ' (3.66) 

t - 1, . « m 

A computer program for calculating the controllability Grammian by numerically 
integrating (3.65) was developed. The output form the program for the example (3,62) is 
shown below and required 0.43 seconds of cpu on the IBM 370. 

M A Tk I fi 

- 0 . 1 0 (J 0 0 0 e 01 
0*0 

0.1 0 00 oof: 01 


^ A T H I X ri 

O.l'^OOOOF Ol -O.lOOOOUtOi 

-o.400oooe 01 o.^ooooot-: oi 

-O.lObOOOE 0^ -0.?b0000F 0<? 

M A I X w 

0.31bl??6F i)2 -0,44450i^E 0^ a.3847S9t: 02 

-O.444S07E 02 0.63617/F 02 -0.444507t; 02 

0.3%47F)9t 02 -0.444b0 7F 02 0.1i9S4bt. 04 

The computer solution of W agrees with the solution (3.63) calculated by hand. 


0.20U000F: 01 u.o 

-O,3U0000E 01 0.0 

0.200000F: 01 -0.2UOOOOt 01 
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Recursive Algorithm 

Suppose the objective is to compute the controllability Grammian for t = T, 21, 

3T, . ♦ NT . Let W(n) denote solution of (3.65) for t + nT and define 

A recursive algorithm for computing W(n), n = 2, 3, . , N from W(l) and O is 
developed below. From (3.64) 

nT+T hT+T nT 

W(n+1) = J Fdt = J Fdt+J Fdt (3.67) 

0 nT 0 

Let 4 = t - nT and from 

4>(-t) = 4>(-«-nT) = 0"4>(-C) 
it can be shown that 

nT+T T n n' 

J F(t)dt = > J F(4+nt)d4 = n ; F(|)d$0 (3.68) 

nT 0 0 

Substituting (3.68) into ,(3.67) and using the definition (3.43) results in 

W(n+1) = n"w(l)ri’' +W(n) (3.69) 

From (3.69) it can be shown by repeated substitution that 

W(n+1) = nV(l)o"' + . . . + OW(l)TV + W(l) (3.70) 

From (3.70) it can be readily proven that 

W(n+1) = OW(n)n' + W(l) (3.71) 

Formula (3.71) can be used to reduce the amount of numerical integration. To compute W(t) ot 
t + NT instead of numerically integrating (3.65) from 0 to NT , only integrate (3.65) and 

- = A<l» , ^(0) - I 

from 0 to T and then use (3.71) 
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3,3,2 Index of Control lab i 11 ty 

When the general vehicle dyrKimics are nonlinear/ then the linear equations (2,9) 
for the dynamic response about the trim solution 6^ are a function of 6^ , Hence, the 
controllability of the linear system (2.9) varies with the choice of the trim solution. 
Quantification of controllability provides a measure for determining the trim solution that 
results in the most controllable linear system. In the previous section, the controllability 
Grammain W(t) at t=T is used evaluate the integral (3.44) for the scalar E which may 
be viewed as the energy expended by the control effectors in returning the vehicle to trim 
during a time span of T seconds. One possible means of quantification is the use of E 
to indicate the degree of controllability. In this section another means of quantification 
is developed. An index of controllability is defined as the ratio of maximum to minimum 
eigenvalues of W(T) or some other controllability matrix. 

The time-invariant linear system (3.39) is said to be controllable , if it is possible to 
find an input u which reduces an arbitrary initial state to zero in finite time T . A 
necessary and sufficient condition for the system to be controllable is that the controllability 
Grammain W(t) defined by (3.43) be nonsingular for some finite t . If W s W(T) is 
nonsingular, then (3,42) defines one of the many possible inputs u that satisfy the definition 
of controllability. Another matrix often used to study controllability is 

t 

P(t) 0(r)BB' |>(r)dt (3.72) 

where 4>(t) is the transition matrix (3.41). The matrix P(t) is related to W(t) by 

P(t) = 4>(t)W(t)iT{t) (3.73) 

and can be identified as the covariance matrix of the state x(t) when u(t) is white noise 
having a spectral density of unity. It follows from (3.73) that the system is controllable 
if and only if P(t) is nonsingular for some finite t . If the system (3.39) is stable, then 
the integral for P(t) exists as t “ and the asymptotic value 

P = lim P(t) (3.74) 

t -* CO 

is the solution to the algebraic equation 
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(3.75) 


AP + PA' + BB' = 0 

It is well known that P(t) or W(t) is nonsingular if and only if the matrix 

K = [B,AB, . . . , a'^'^’b ] (3.76) 

has rank k = order of the system. The rank of K is equal to the rank of kxk symmetric 
matrix 

Q = kK' = BB' +ABB'A' + . . . + A*^‘^^BB'(A')*^”^ (3.77) 

which is more convenient than K for testing controllability. 

Indices of Controllability 

The necessary and sufficient conditions for controllability of a time— invariant system is 
that a certain matrix be nonsingulor. Possible choices of the test matrix that ore symmetric, 
positive-semidefinite include W(t) , P , and Q . This controllability is a property that a 
given system theoretically either possesses or does not possess. In practical applications, 
however, there may be instances in which a system may be nearly uncontrollable in the sense 
that certain initial states moy be much harder to reduce to zero than others. Evidence of such 
situations is that the matrices tested for controllability are nearly singular, i.e., poorly- 
conditioned. It is thus appropriate to use conditioning of a relevant matrix as an index of 
controllability. A useful measure [4, 5 ] of the conditioning of a matrix F is 

MF) = 11 F 11 • II F‘’|| ' (3.78) 

where jj F 1| denotes the norm of the matrix F defined by 
11 F II = sup 11 Fx II ^ 

l|x||=1 

where || x j| is a suitable vector norm. When the Euclidian norm, i.e. 

• II X II = 

is used, then, fora symmetric matrix F, 
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I X / X . I 7Q\ 

' max min ' 

where and \ are me eigenvalues of largest and smallest magnitude, respectively. 

Clearly k(F) s 1 and reaches the lower limit only when I X I = I X I . i.e. 

' ' max ' I min ' 

when oil eigenvalues ore equal in magnitude. The condition k(F) ^ 1 also holds for other 
norms, as shown in [5^. 

The quantification of controllability (and/or observability) was considered earlier 
by several investigators. Kalman, Ho and Narendra [5] considered using the trace or 
the determinant of the inverse of the controllability matrix as indices of controllability, 
and Johnson [^7] considered the determinant as an index of controllability in greater detail. 

The shortcoming of the earlier indices of controllability is that they depend on 
the scale of the variables used in the problem. For example, multiplying earch control 
variable by a constant say c , is equivalent to multiplying the B matrix by the same 

constant and hence the controllability matrix Q as defined by (3.77) or P as defined 

2 -1 -1 ~2k 

by (3.74) is multiplied by c . Hence the trace of P or Q is multiplied by c 

On the other hand the conditioning number is obviously independent of a scale change, 
either of the control variables or of the state variables. The conditioning number, how- 
ever, does depend on the choice of stote variables, as the following examples Indicate, 

Example - Consider the system having the transfer function 


H(s) = 


Y(s) _ s + a 

(s + l)(s + 2) 


It is clear that if a = T or 2 the system is either not observable or not controllable or 
both. The objective of this example is to show the behavior of the controllability 
index as a -* 1 or 2 . 


In order to examine the controllability and observability of the system it is 
necessary to define a suitable set of state variables. In this example the state variables 
are defined as those of two canonical forms. The Jordan normal form and the companion 
form. 


40 



Jordan Form - The Jordan form con be obtained by expanding H(s) in partial fractions: 


H(s) - 


a - 1 
s + 1 


g - 2 
s + 2 


Two block diagram representations of H(s) are given in Figure 3.5. For Figure 3.5(a), 
the state and output equations are 


X = -X + u 


-1 

0 

rn 

1 1 

A = 



B = 1, 



0 

-2 

N 


X = -2x + u 
2 2 


y = (a-l)x - (a-2)x C - [Q"1/ “(o“2)] 

12 

For Figure 3.5(b) the state and output equations are 


X = -X + (a-l)u 


-1 

0 ' 


a-1 ’ 

1 1 

A = 

. 0 

-2 

B = 

- -(a-2) 


X = -2x - (a-2)u 
2 2 


y = X +x C = [ 1 1 ] 

12 

It is thus seen that the A matrix of both representations are identical and also 

. Moreover, B' of Figure 3.5(a) equals C of Figure 3.5(b) and C of 
Figure 3.5(a) equals B of Figure 3.5(b). Hence it follows that observability of Figure 3.5(a) 
corresponds to controllability of Figure 3.5(b), and vice“versa. Accordingly, examining 
the controllability of Figure 3.5(b) is equivalent to examining the observability of 
Figure 3.5(a). 


The controllability matrix K of the system of Figure 3.5(b) is 


K = 


a-1 

-(0-2) 


-{o-D 

2(a-2) 
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Hence 


Q = KK' = 


2(a-l)^ 

-3(o-1)(a-2) 


-3(a-l)(a-2) 

5(a-2)2 


The characteristic equation of Q is 


- X(trA)+ lA] = 0 


where trA = 2(a-l)^ + 5(a-2)^ 

1A|= (a-l)V2)^ 


There is o characteristic root at \= 0, for a - 1 or a - 2, and these are the 

values of a for which the system is not controllable, as expected. The condition 

number of Q / as defined above, is 


k(Q) = 


trQ + v(trQf - 4 JQ 


trQ - y<trQ)2-4| Q 


A curve showing the behovior of k(Q) vs the parameter a is shown in Figure (3.6). It is 
observed that k(Q) tends to infinity as a - 1 or as a - 2. It is interesting to note, 
however, that k(Q) reaches (local) minima of 37.9 at a =1.61 and a =3.72. This would 
suggest that if a were adjustable, the controllability (or observability) can be optimized, 
in the sense of minimizing k(Q) by using o=1.61 or o = 3,72. 


Instead of k(Q) we can determine k(P) after solving for P by use of (3.77) 
The solution of the latter is 


P = 


(a-l)2 

(a-l)(a-2) 


(a-.l)(a-2) 

( 0 - 2)2 
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-2 


(a) (b) 


FIGURE 3.5: JORDAN CANONICAL FORMS OF 

TRANSFER FUNCTION IN EXAMPLE 



(°) (b) ' 


FIGURE 3.8: COMPANION FORM OF TRANSFER FUNCTIONS 
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whence trP = 


( 0 - 1)2 („_ 2)2 

|P|= ^(o-l)2(a-2)2 

The resulting curve for k(P) is also shown in Fig« 3.6. It is observed that k(P) attains 
minima of about 34,0 at aj5^-i,5 and a j5-fl,4. 


It is noted that the minimum value of the conditioning number Is almost equal for 
P and Q and one minimum occurs (as expected) between a = 1 and a = 2. The 
locations of the other minima are quite different, but the general shapes of the curves are 
remarkably similar. 


It is of interest to examine the effect of adding another independent input on the 
controllability of the system. Suppose, for example, another input say u^ was added to 
the first state, resulting in the equations 

X| = -xj + (a-l)u] + U 2 

X 2 = - 2 x 2 " (a-2)u] 

The corresponding B matrix Is now 


B = 


a-1 


L-(a-2) 


1 

0 


The controllability matrix is now 


K = 


a-1 

L-(a-2) 


-1 

0 


-(a-1) -1 

2 (a -2) 0. 


and 



Q = KK' = 


2(a-l)2+2 

-3(a-I)(a-2) 


-3(a-l)(a-2)‘ 
5(0 -2)2 
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likewise 


P = 


(o-1)2 + 1 


(o-I)(a-2) 


-(a-l)(a-2) 

5— 

( 0 - 2)2 

"T“ 


P and Q are now singular For only one. value of a , namely a = 2; obviously 
X 2 is no}- controllable for a = 2, 


The curves of k(P) and k(Q) are shown in Fig. 3*7. It is noted that the addition 
of input U 2 has the effect of reducing the conditioning number for all values of a , as 
would be expected. 

Componion Form -- Two alternate companion forms that realize the transfer function H(s) 
are shown in Fig, 3,8 (a) and (b). The corresponding matrices are as follows 


Figure 3, 8 (a) 

A = 


Figure 3.8(b) 

A = 



0 

-2 



B = 


1 

a -3 


B = 




c= [' o] 

’] 


Since the C matrix of Fig. 3.8(a) is independent of a it is natural to examine the behavior 
of this realization for controllability. Likewise, it is natural to examine the realization of 
Fig ,3,8 (b) for observability. 

For the system of Fig. 3.8 (a) it is found that 


1 a - 3 

a — 3 —3a + 7 


whence 



" a^- 6a + 1 0 

-3a2 + 17a-24 


-3o2 + 17a-24 
10a^-48a58 ^ 
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and hence tr Q = lla^-54a+68 

|Q| = (a-1)2(a-2)2 

Solution of (3.75) for P gives 


' a ^+2 

nr" 

I 

_ "2 


2 

2 


a^-da + ll 

"n 


with 


tr P = ^(a^-4a+8) 
jPl =^(a-l)2(a-2)» 


Curves Showing k(P) and k(Q) os functions of a are given in Figure 3.9. |t ij 
noted that although local minima occur for both k(P) and k (Q) for 1 sa ^ 2 , the minima 
attained exceed 1000 and hence would indicate that operation with a in this interval is 
undesirable. A very sharp local minimum in k(Q) of about 1.4 occurs at a , and 

would indicate that operation at this value of a is, in a sense, optimum; k(P) on the 
other hand does not have any other minimum, but tends to unity as a -* ± <», This corresponds 
to the case in which the "feedforward” gain (to ) is negligible in comparison to the 
direct gain (a -3). 
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3.4 OPTIMUM CONTROL APPROACH 


3.4.1 Optimum Control Computation 

If the general control problem described in Section 2. 1 can be formulated as an optimum 
stochastic control problem for a linear process with a quadratic performance criterion than a 
linear feedback system can be designed to solve both the trim problem and the dynamic response 
problem. The theory to compute such a feedback system is developed In this section and will 
be applied in Section 4.3 to the lateral control of the Space Shuttle. 

The linear stochastic optimum control problem with bias inputs ‘is defined by the following 
equations in vector-matrix notation 


Process Dynamics: 

X = Ax + Bu + Cz + v 
Efv} = 0 

Observation Equation: 

y = Hx + w 
E{w} = 0 

Performance Criterion: 


2 = constant 
E{vv'}=V 


E{ww'] = W 


9 

J(u) = E{ J(x'Qx +a o'.Ru)ds l y(7) for T^t] 
t 


‘ O' = scalar parameter 


where 


X = state vector 
u = control vector 
y == output vector 
z = bias vector 

V = input noise vector to process dynamics 
w = sensor noise vector 


(3.80) 


(3.81) 


(3.82) 


Equation (3.80) is identical to (2, 10) except the vector of deflection angles is denoted by u 
instead of 6 , The stochastic optimum control solution is denoted by u in order to distinguish 
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It from the control solution 6 obtained by solving the trim problem. This distinction is 
helpful in the next section when the correlation between u and 6 is developed. In the 
usual problem formulation, the scalar parameter a is not present in (3.82) since it can be 
incorporated in the R matrix. In this case, however, the scalar parameter a is useful 
in deriving the correlation between u and 6 . 


If the "bias term" Cz was not present, then the optimum control problem defined by 
(3.80) - (3.82) would be in the standard form. By defining z as part of the state vector, 
(3.80) - (3,82) may be rewritten in the standard form. The resulting augmented dynamics are 


X = ^x + Bu+v 

(3.83) 

y = H X + w 

(3.84) 

J = E J(x' Q X + u'Ru)ds I y(r) for t ^ t] 

(3.85) 


t 


where 



X 


A 

C 


B 


Q 

1 

0 

X = 

z 

j : = 

0 

E 

B = 


5’ = 

.0 

0 ^ 












V 


" V 

o" 

/• 


1 



V = 

f 

V = 

0 

* 

3*_ 

H 

— 

H Oj 




The solution to the optimum control problem defined by (3.83) - (3.85) is given by the 
equations 

Deterministic Quadratic Optimum Control ; 

u(t) = - Fx(t) (3.86) 
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where 


(3.87) 


F = ( 1/aVVM 

MA" + 7f'M-0/o^)MBi?''B'M+^ = 0 (3.88) 

Kalman Filter: 

■ ■ ■■ ^ -A A 

X = Ax + Bu + K(y-Hx) (3.89) 

where 

K = pTT'W"’ (3.90) 

0 = yp + p/ST' - pff'w"’HP +7 (3.91) 

The feedback control system defined by (3.86) - (3.91) is divided into two parts in tandem. 

A 

First, a Kalman filter computes the optimum estimate of the augmented state x from the 

A 

sensor measurements y . Next, feedback gains multiply the estimated state x to yield the 
control signal. In the event that the augmented state vector x can be measured perfectly, 
i.e., 

yHx 

then, the Kalman filter is not required. In this cose the control system is defined by (3.86) - 

A 

(3.88) where x = x . 

Partitioning theaugmented state vector into x and z simplifies the equations (3.86) - 
(3.91) for the control design. The deterministic quadratic optimum control is considered first. 


By partitioning the matrix M according to 


M = 


M, 

iMj M3 


the optimum control solution (3.86) can be rewritten ds 


u(t) - u (t) + u (t) 

X z 

where 

u (t) = - l/g2R‘’B'M,x(t) = - F 5(t) 

X I X 

u^(0 = - = - F z(t) 


(3.93) 
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The symmetric matrix is the positive definite solution of 


M 


-hr 


+A'M^ - l/a^M^BR B'M^ + Q = 0 


(3.94) 


and the matrix is computed from according to 


= - (A'-l/a2M|BR‘’B')"'M C 


(3.95) 


In the derivaHon of (3.94) and (3.95) it is assumed that E = 0 in S' which corresponds 
to the assumption z = constant . 

Similarly, by partitioning the matrix P according to 


P = 


P P 

1 2 


P' P 

3 


the equations (3,89) - (3.91) for the Kalmon filter become 


X = Ax + Bu + Cz + K^(y-Hx) 


(3.96) 


and 


z = Ez + K (y~Hx) 


K = P.H'W 
X 1 


k = PIH'W 
z 2 


-1 


-1 


(3.97) 


he partitioning of the P matrix does not simplify the computation of the submatrices P^ and 
' as in case of the matrix M . Hence, P- and P 2 are computed by solving (3.91) for the 
KJsitive definite covariance matrix P . In the computation of P it is assumed E 0 and 
0 . If 0 then ^2 “ ^3 ~ ^ that the bias disturbances z can be 

letermined perfectly which is not realistic- A small amount of damping (E ^ 0) is included 
n the noise model of bias disturbances in order to yield a finite value of P^ . 
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3.4,2 Correlation Between Trim Solution and OpHmum Control Solution 

There is a relationship between the optimum control approach and the trim, control 
approach. This relationship relates the optimum steady state control value u{«) to the trim 
solution 6 for the case when the control weighting matrix R in the performance criterion 
(3,82) of the optimum control approach and in the performance criterion (3.19) of the trim 
control approach are the same. 

The derivation given below is for the case of complete state feedback for which (3.92) 
holds. It appears that the proof extends to the more general case in which the optimum 
control system includes the Kalman filter to estimate the state. A detailed proof, however, 
has not been developed for the more general case. 

Substituting (3.93) into (3.80) yields for the case of complete state feedback the closed 
loop dynamics 

x= Ax +Cz (3.98) 

where 

A = A - 

C = C - 1/(7*BR‘’b'M, 

Since the matrix A is asymptotically stable, setting x = 0 in (3.98) results in the formula 

x(“)=-A"’cz (3.99) 

for computing the steady state value of the state vector. In turn, substituting (3.99) and 
(3.95) into (3.93) gives that the steady state value of the control vector is 

u(“ ) = u (“) + u (“) 

where 

u (») = R"'b'M, (a2A)"’cz 

u^(») = R‘’b'((t2a')‘'m,Cz 


(3.100) 

(3.101) 

(3.102) 
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Next we will consider how u(») varies with the scalar parameter a >n the 
performance criterion (3.82). In particular what is the limiting value of u(o°) as a 
approaches zero. In determining the limiting solution, we must take into account the variation 
of the matrices and with a • A solution of (3.94) is sought in the form of a series in 
ascending power of o* J 


= Nq + (tN^ +aN2 + . . 


(3.103) 


In papers by Friedland [6] and Hutton [9], it is shown that the following equations: 


NoB=0 


(3.104) 


NqA + A'N. + Q - NjBR"’b'Nj = 0 


N|A+A'N, - N2BR"'b'N^ -NjBR'^B'Nj = 0 


(3.105) 


must be satisfied if (3. 103) is a solution to (3.94). The above equations ore formed by 
substituting (3. 103) into (3.94) and equating matrix coefficients of like powers of cj . By 
matrix manipulations of (3. 104) and (3. 105), it is shown in [B] that Nq is the positive semi- 
definite solution of 



N AU - B(B'QB)'’b'Q] +Cl - QB(B'QB)‘'b']A'N|j + Q - QB(B'QB)"’ B'Q 
- NqAB(B’QB)"'b'A'N^ 


(3.106) 


After solving (3.106) for , we can solve (3.105) for the positive semi-definite matrix 


Consider the asymptotic value of 

(o^A)"’ = (ct^A - BR”’b'M )■' 


(3.107) 
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as a approaches zero. For a I ! nonzero 0 , the matrix is positive definite. From 

(3. 104), the matrix M is positive semi-definite at 0=0. However, if the first term 
2 ; 

0 A in (3. 107) decays to zero more rapidly than the second term, then 


(o^A)"’ - - osa^ 0 


(3.108) 


provided BR is positive definite. Substituting (3.103) into (3. 107) and using (3.104) 
gives that 


C(7^A]“’ = [a^A - BR‘'b'(otN^ +o^N2 + , . . )]"’ 


(3.109) 


The dominant I’erm in (3. 109) is BR B'N. which is derived from the second term in (3. 107) 
and indicates that (3. 108) is valid. 


Substituting (3. 108) into (3.95) and (3.98) gives 


lim /VL = - (BR'^B')“^C 
a-+0 ^ 


(3.110) 


lim C — 0 
a -^0 


(3.111) 


Further substituting (3.108) into (3.101) and (3.102) and using (3.111) yields the results 
lim u(“)= lim u (“ ) = - R“^B'(BR"^B')’‘^Cz 


0-»O 


a -»0 


lim u (“) = 0 

V ' * 

CT -*0 


(3.112) 

(3.113) 


The trim control problem is to find the set of controls 6 satisfying 
0 = Cz + B6 


(3.114) 


and minimizing the performance index 
J = 1/2 e'R6 


(3.115) 
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The solution to (3. 114) and (3.115) is 

6 = -r’B(Br'B'r'cz (3.„<5) 

Comparing (3.1 16) to (3. 1 12) provides the fundamental result that 

6 = limu(~) if R = k^R (3.117) 

a -»0 

where k is an arbitrary scalar. Thus the steady state value of the optimum control solution 
in the case of unlimited control authority (control weighting matrix R in the performance 
criterion goes to zero) is egual to the trim solution provided the relative control weighting 
matrices are the same in both cases. This provides a correlation between the optimum 
control solution and the trim solution. 
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4. SPACE SHUTTLE CONTROL 

Control of the Space Shuttle is studied during ascent when more control effectors 
are available than required. The analytical methods developed in Section 3 are applied to 
the lateral control problem. An illustration of the Space Shuttle configuration, given in 
Figure 4, 1, shows the two aerodynamic surfaces and five rocket engines available for control. 
For purposes of later reference these controls are identified as follows: 

1) top orbiter rocket engine 

2) right orbiter rocket engine 

3) left orbiter rocket engine 

4) right solid rocket motor 

5) left solid rocket motor 

6) aileron, 

7) rudder 

By varying the angular position of these^controls, seven independent means of lateral control 
are achieved. But only three independent controls are required, leaving four redundant 
controls. If the solid rocket motors (SRM) are not gimballed then the number as independent 
controls is reduced to five, leaving two redundant controls. The results in this report are 
for the latter cqse. However, the equations and computer programs used to perform the 
calculation of the control deflections include the possibility of gimballing the SRM. 

4. 1 Space Shuttle Dynamics 

A mathematical model describing the lateral motion of the Space Shuttle is given in 
this section. This description entails an extensive number of the parameters defined in 
Appendix B tpgether with a tabulation of their numerical values . 

The set^of. differentia I equations describing the translational and rotational motion of the 
vehicle are based on summing the forces and moments along the body axes of the vehicle *, 

The body axes are defined as a Cartesian coordinate system fixed to the vehicle and whose 
origin is located at the center of mass as shown in Figure 4.2. The attitude and rotational rate 

* The notation and definitions used for the aerodynamic terms in the report are in accordance 
withC 1], 
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of the vehicle are defined by the Euler angles and the components of the angular velocity 
vector along the three body axes. Specifically 

(p = roll angle 
0 = pitch angle 
ij) = yaw angle 
p = roll rate 
q = pitch rate 
r = yaw rate 

1 

Figure 4.2 Body Axes and Notation 




The positive directions are as shown in Figure 4.2. Upper case letters denote the total motion 
(Nominal and Perturbation) of the vehicle. The velocities, forces, and moments about the 
three body axes are defined as follows; 

U, X = forward velocity and force 

V, Y = side velocity and force 

W, Z = downward velocity and force 
L = rolling moment 

M = pitching moment 
N = yowing moment 


The kinematic and dynamic equations describing the lateral motion of the vechile are 
Y = m[V +RU - PW - geos e sin$] 


L=JP-I R+QR(j-r)-I PQ 

X xz z y xz 


(4.1) 


N = -J P + J R + PQ(J -J ) + J QR 

xz z y X xz 


and 


$ = P + Q sin ^tan 0 + R cosOtan© 


(Qsin^+Rcos$)sec6 
where the moments of inertia are defined by 

^x " J* 

I = J(x^+y^)dm 


- J (x^-tz^) dm 

I = r xzdm 
xz 


(4.2) 


(4.3) 


For this investigation no data was available on I , thus the approximation I — 0 is 

X Z , X z 

used. The (total) vehicle motion modeled by (4. 1) and (4.2) can be partitioned into nominal 
plus perturbation rrwtion by substituting 


U = U +u 
o 


P = P +p 
o 


© = © + 0 

o 


V = V +v 
o 


Q = Q + q 
o 


$ +<n 
o ^ 


(4.4) 


w = W + w 
o 


R = R +r 
o 


^ = 'I' + i/) 
o ^ 
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where the capital letters with subscript “o" denote the nominal motion and the lower 
case letters denote the perturbation motion. For the nominal motion along the trajectory 
it is assumed that 

U?^0 P=0 ey^o 

o o o 

V=0 R=0 ^=0 (4.5) 

o o o 

W 7^ 0 Q 0 ^ = 0 

o o O 

The nonzero values are tabulated in Appendix B for each of the twelve flight times along 
the ascent trajectory for which the perturbation motion is to be studied. Substituting (4.4) 
and (4.5) into (4. 1) and (4.2) results in the following linearized equations of motion for 
small perturbations from the nominal trajectory: 

Y = mCv + U r-g cos © ml 
o o^ 

I = I p +(J -I )Q r 
X 2 y o 

N = J r + (J-J )Q p (4.6) 

z y X o*^ ' ' 

O = p + Q tan © <a+ tan © r 
^ o o^ o 


J) = (Q o + r) sec © 

Adding the equation y = v to (4.6) and reqriting in the state space formulation 
results in the vector-matrix equation 

X = Xx + Bf (4.7) • 

where the state and forcing vectors are 

X = Cy / <p / / V , p , r ]' 

and 

f = [Y, L , N]' 

The constant matrix B has the form 




where 


A = Diag [ m , J , J ] (4. 

The forcing vector f represents the lateral forces and moments acting on the vehicle and 
can be modeled by 

TyI Fy y yIFvI 

V p 


f = U = L 


p +B6 + Ci 


n ^p ^JL/J / 

aerodynamic forces and control bias disturbance 

moments forces forces and moments 

and 

moments 

Substituting (4.9) into (4.7) gives the desired vector-matrix equation for the dynamics of 
the vehicle 

X = Ax + Ba + Cz (4. 


where 


0 0 0 1 


°22 

0 

0 

1 

°26 

°32 

0 

0 

0 

°36 

®42 

0 

^44 

°45 

“46 

0 

0 

^54 

®55 

“56 

0 

0 

^64 

°65 

“66 

= Q tan© 





a_, = tan© 

40 O 


= Q sec© 
32 o o 


a«. = sec© 
36 o 


°42 " 


'44 

= Y 

V 

.45 p 

°46 

II 

1 

c 

o 

'54 

= L 

V 

a__ = L 
55 p 

°56 

= L +Q (J 
r o y 

'64 

> 

Z 

II 

o.c = N +Q (J -I )/J 
65 p o' X y'^ z 

°66 

- N 

r 
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and 


3 



r 


B = 

0 


-1^ 



A B 

' 


3 

3 



(4J2) 


In the remainder of this section, the formulas for computing the matrix elements in (4.9) 


which are required for (4. 10), are developed. 

The formulas for the matrix elements corresponding to the aerodynamic forces and 
moments ore 


Y 

= QC VU 

II 

o 

Y 

= QcC /2U 

V 

0 

1 

p 

r 

yr o 

L 

= Q C ,U 

L = Q bC, /2U 

L 

= Q bC, /2U 

V 

1 

X 

1 

0 

p X ip o 

r 

X tr o 

N 

V 

where 

= Q C VU 
z n0 o 

N = Q bC /2U 

p z rp' o 

N 

r 

= Q bC /2U 

z nr o 

Q 

~ qS/m 

q = 

= qSb/I 
dynamic pressure 

Q 

z 

= qSb/j^ 


U = nominal velocity in x-direction 
o 

S = reference area 

b = reference length 

c = length of mean aerodynamic cord 


(4.13) 


Next the expressions for the forces and moments generated bygimballing the rocket 
engines are derived. The location and nominal direction of each rocket engine with respect 
to the Cartesian coordianate system fixed to the vehicle is shown in Figure 4.3 * . The rocket 
engines are numbered 1 through 5 as indicated in Figure 4,3 and in agreement with the list 
of controls at the beginning of Section 4. Let x^, y^ , denote the coordinates of the 


^ The location of SRM was not included in the information received from MSFC. This data 
was not required since it was assumed the SRM could not be gimballed. However, the equations 
and corresponding computer programs include the posibility of gimballing the SRM. 
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vehicle center of gravity where y = 0 , The (position) vector from the center of gravity 

C9 

to the tth rocket engine is, therefore, 

[x +x , y. , z +z 1 (4‘. M) 

eg t ' t eg t 

The thrust vector with magnitude F has the components 

Tjf j 

(forward) . ~ 6^cos cos 9^ sin " sin 0^os 

(sideward) = F^cos 6^sin +cos6^cos j/^6^^^-sin 6^sin0^^^^) (4-15) 

(downward) Z = F (sin 6. + cos 6.6 J 

1/ L> if ep^ 

where the angles defining the direction of the thrust vector are 

9^ = nominal pitch angle of the tth rocket engine. 

t ■ 

ll) = nominal yaw angle of the tth rocket engine. 


6 . = pitch deflection of the tth rocket engine. 

epi ^ 

5 = yaw deflection of the ttb rocket engine. 

eyt 

as shown in Figure 4.4. The arrows in Figure 4.4 indicate the directions of positive angles. 
The nominal directions of the rocket engines are shown in Figure 4.3 and listed in Table 4, T, 
The derivation of (4. 15) assumes that the deflection angles are small . 

Table 4. 1 Nominal Directions of the Rocket Engines 


index i 


1 

- 18“ 

0 

2 

- 12® 

- 3.5® 

3 

- 12® 

3.5° 

4 

0 

- 15® 

5 

0 

15® 
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Figure 4.4 Angular Direction of Thrust Axis 



z 
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The moments induced by the ^-th rocket engine are given by the cross product of the petition 
vector (4.14) with the thrust vector which results in 


(roll) =y^Z^-(z^g+z^)Y^ 

(pitch) 

(yaw) Nj = (x^g+x^)Y^-y^X^ 


(4. 16) 


Substituting (4.15) into (4. 16) expresses the moments as a linear function of the deflection 
angles. 

Having derived the general equations for modeling the rocket engines, the next step 
is to derive the equations corresponding to the term B6 in (4.9). 

The elements of the control vector are 

6 . = 6 


1 


eyl 


1 


®ey2 2'^®ey3'^®ey2^ 
®3 " ®ep3 " 2^®ep3‘ ®ep2^ 
= ®ey4 ^ 2^®eyS'^®ey4^ 
®5 " ®ep5 ^ ?^®ep5‘®ep4^ 


6 , = 6 

6 a 


(4.17) 


«7 


= 6 


where the deflection angles are defined as follows: 

= yaw angle of top orbiter engine 

angle of right orbiter engine 

6 « = pitch angle of right orbiter engine 

ep4 

= yaw angle of left orbiter engine 
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6 - = pitch angle of left orbiter engine 

epo 

^ey4 ^ angle of right SRM 

^ep4 ” angle of right SRM 

6 r = yaw ongle of left SRM 
eyo ' 

6 - = pitch angle of left SRM 

ep5 


The elements of the constant 7x3 matrix 



hi 

bi 2 

h 3 

h 4 

*^15 

h 6 


II 

hi 

h 2 

hs 

■=24 

hs 

*’26 

^27 


hi 

h 2 

hs 

h 4 

^35 

*=36 

^37 


are computed from the following set of formulas; 
b ^ ^ = F cos 18® 

b2, = -F(z,-Zcg)cosl8“ 

b«- = F(x,-x ) cos 18® 

31 1 eg 

b ^2 = 2F cos 12® cos 3. 5° 

b«« = -2F(z -z ) cos 12® cos 3. 5° 

b__ =2F[(x_-x )cos3.5‘’ -y„sin3.5° ]cos 12® 

32 2 eg 2 

b., = 2F sin 12° cos 3.5° 

^23 = 2FCy2Cos12°-(z2-z )sin12°sin3.5°] 
bgg = 2FCy2cos3.5°+(x2-x )sin3.5° ]sin 12° 


{ 4 . 18 ) 


( 4 . 19 ) 
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‘’14 = 2Fsrm~' 

‘’24 

‘’34 = 15» -y^sin 15“ ] 


‘’15 = ° 


'25 


= 2F 


SRM^4 


‘’35 = ° 


‘’16 = 




^26 = 


36 


qSb ,(C,. ) 

^ ref .to eg 
a 

qSb iC - ) 

^ ref n6 eg 
a 


(4. T9 continued) 


‘’17 ‘’^^6 

' r 

‘’27 = ^cg 

r 

S7 = >cg 

r 


^*-46 ‘eg Se ^^cg"*mrp‘/‘ref 

a a ' a 

(C . ) = G . - C ^ (x -X )/b , 

no eg no yO eg tnrp ref 
a a ' o 

(4.20) 

^^ 06 ^‘cg = "^Ss/^cg ■ *mrp‘^‘’ref 

^^'nB^^eg " ^n6^ " ^y6^\g~^tmp^ ‘^Ket 

The formulas in (4.19) are grouped by column. The tth column of the B matrix in (4.9) 
defines the values of Y , L., N corresponding to 6^ « The formulas for the first five 
columiis are derived from (4.14) - (4.17) . The last two columns corresponding to the aileron 
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and rudder, respectively, are computed using the standard formulas for aerodynamic control 

surfaces. The data for the stability derivatives received from MSFC were with respect to the 

moment reference point located at x , y , z where y = 0 . The translation 

mrp mrp mrp mrp 

of data from the moment reference point to the center of gravity is given by (4.20). 

The force and moments in (4.9) due to the bias disturbances is modeled by the term 
?z . The elements of the vector z or bias inputs are 

= ^ = side slip angle due to a steady side wind 

z = T = roll bias torque due to SRM misalignment 
^ *'b 

= T = yaw bias torque due to SRM misalignment 
The constant 3x3 matrix C" has the form 




Si 1 0 


Si ° ' 


where the elements in the first column are computed from 


n 


= q 


sv 


Si = 

Si = '<sb(c;;^cg 
S/3 ” Ss ^ ^^y/3 

S/3 “ S/3 ^^ne^FORWARD 


(4.21) 


(4.22) 
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The last three equations in (4.22)account for the change in the stability derivatives due 
to the addition of a pair of dorsal fins to the Spoce Shuttle configuration os indicated in 
the sketch below. 



To summarize, the lateral dynamics of the space shuttle is governed by the vector 
motrix equation (4.10), The coefficient matrices A , B , C in (4.10) are computed using 
(4.8), (4.9), (4. 11) - (4.13), (4.18)- (4.22). The values of the parameters required by 
these equations are given in Appendix B, 
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4,2 TRIM PROBLEM AND SOLUTION 


When bias disturbances generate forces and moments that cause the vehicle to deviate 
from the nominal ta jectory, the rocket and aerodynamic controls must be deflected in such 
a way as to counterbalance these forces and moments. For the lateral trim problem of the spoce 
shuttle, the bias disturbances are due primarily to steady side winds and SRM misalignments. 

The state vector x in (4.7) defines the deviation of the vehicle from the nominal trajectory 
in the lateral -direction. Hence, the trim condition is to maintain x = 0 . On substituting 
X = 0 into (4.10) one finds that the trim solution ^ must satisfy the matrix linear equation 

0 = B6 + Cz (4.23) 

For a given value of the bias vector z , (4.23) represents six equations in seven unknowns. 
However, the equations are not all linearly independent. From (4.12) the first three equations 
are identically zero independent of 6 and the last three equations have the form 

0 = A^'^a + a’^z (4.24) 

where A is the diagonal matrix defined by (4.8). Premultiply (4.24) by A gives 

0 ='S'6 + Cz . (4.25) 

which is equivalent to setting Y = L = N= 0 in (4.9). In other words, (4.25) states that 
the trim control must provide zero net side force, rolling moment, and yawing moment in the 
presence of a steady side wind and SRM misalignments. Replacing (4.23) by (4.25) has reduced 
the number of trim equations from six to three. In terms of the notation introduced in Section 
2. 1, the dimensions of the trim problem are 

m = 7 : number of controls 

n = 6 : number of state variables 

n = 3 : number of linearly independent trim equations 

In order to determine the optimum trim solution, a performance criterion of the following 
form was selected; 
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r(6)=5 ^ 


r (1-C0S5 ) 

7l ^ 

6 


(4.26) 


where 


max 

q 

S 



Maximum deflection angle allowed for the ^th control because of physical 
limitations or excessive hinge moments. 

Dynamic pressure 

Reference area corresponding to the drag induced by the t^h control 
Coefficient of drag corresponding to the ^th control. 


The numerical values of the above parameters is given in Appendix B. 


iT 

The seven components of the vector 6 of control deflections are defined according 
to (4#17). The first term in (4.26) penalizes the movement of the actuators for trim in order 
to leave maximum flexibility for dynamic response. The second term in (4,26) penalizes the 
thrust loss(gain) caused by gimballing the rocket engines away from their nominal position. 

The third term in (4,26) penalizes the thrust loss due to drog caused by deflecting aerodynamic 
surfaces. 

Substituting the approximation 
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into (4.26), the performance criterion can be written as the quadratic form 

r(6) = j6'R6 (4.27) 

where R is a diagonal matrix whose elements are given by 

R = w2 /6^ +W2 (4.28) 

it \t 2h 

R =wVe^ +WMqS.C^)2 i = 6,7 (4.29) 

it \t trr\ax it i 

The fourteen (relative) weighting factors W and W are selected by the user to achieve 

2i 

the best performance within the restriction imposed by the problem. This best performance is 
a judgement evaluation unless additional criteria are used. 

The lateral trim deflection angles are the solution to the optimization problem defined 
by (4.25) and (4.27). The objective is to solve the trim problem for the maximum expected 
values of sideslip angle and for different combinations of roll and yaw misalignment torques 
that encompass the worse case situation. The sideslip angle is computed from the mean, side 
wind velocity and the vehicle velocity according to 

jS = sin \v^) 

The values of V and V for each of the twelve trajectory points are listed in Appendix B 
and result in the values of sideslip angle listed in Table 4.2. Plotting the values of $ 
as a function of flight time yields the sideslip profile shown in Figure 4.5. Eight different 
combinations of yaw and roll bias torques due to SRM misalignments were provided by MSFC 
for studying the trim problem and these are listed in Table 4.3. 

A computer program entitled TRIMS for computing lateral trim of the Space Shuttle was 
developed. The TRIMS program solves the trim problem given by (4.25) and (4.27) using the 
numerical methods described in Section 3.1 . The program user can select either the steepest 
descent method or the Newton-Raphson method at execution time. Although the trim problem 
given by (4.25) and (4.26) is linear, these numerical methods have the capability to solve the 
nonlinear problem. The TRIMS program is coded to facilitate changes in the trim problem 
including the replacement of the linear trim problem by a nonlinear trim problem. 
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Table 4*2 Sideslip Angle for Different Flight Times 



flight 

tinie 

(sec) 

6 

(rad) 

8 

(deg) 

25 

.02096 

1.201 

40 

.05996 

3.436 

50 

.07887 

4.519 

60 

.09942 

5.696 

65 

.10642 

6.097 

70 

.11124 

6.374 

75 

. 1 1635 

6.667 

80 

.11404 

6.534 

90 

.06169 

3.535 

100 

0 

0 

110 

0 

0 

140 

0 

0 




Table 4.3 Bias Torques Caused by SRM Misalignment 
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Ftgure 4.5 Sideslip Angle v$ Firght Time Due to Mean Wind Disturbance 



The formulas developed In the previous section for computing the matrix elements in 
(4.25) are coded into the TRIMS program. The numerical data required by these formulas 
and tobulated in Appendix B is also stored internally in the program. Similarly the formulas 
(4.28) and (4.29) used to compute the performance criterion (4.27) are coded into the 
program together with the required numerical data. Only those input parameters with values 
that are likely to vary from run to run are entered as input data at execution time. These are 
the fourteen weighting factors w^^ and w^^ in the performance criterion and the values of 
the roll and yaw bias torques, and / respectively, A more detailed description of 
the TRIMS program including flowcharts, listing, instructions showing how to use the program 
is given in Appendix C, 
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The trim angles for the eight different combinations of yaw and roll SRM bias torques in 
Table 4,3 were computed in a single run of the TRIMS program. Each case entailed computing 
the trim ongles for the twelve trajectory points or flight times which totals to 96 trim solutions. 
The total cpu time was 5.29 seconds on the IBM 370/165 computer which averages to 0.055 
second per trim solution. For this run the second order gradient method and the weighting 
factors in the performance criterion were chosen to be 

^ ^ 1 3000 for t/ 6 
[4000 for i=6 

^ t = 1 , . . . . 7 

The laterol trim solutions for the eight cases in Table 4.3 are shown in Figure 4,6 
where the trim angles are in degrees. The trim angles 6^ and 6^ for the SRM engines are 
always zero since in the current TRIMS program the SRM engines are not gimballed. However, 
the provision for gimballing the SRM engines has been included in the development of the 
development of the TRIMS program. 

For most of the trajectory points in Figure 4.6, especially those with high dynamic 
pressure, some of the deflection angles exceed the allowable limits by an order of magnitude. 
This indicates that the Space Shuttle configuration does not have sufficient control authority 
to meet.the trim conditions L = M = 0 when the SRM engines are not gimballed. 

A check of the TRIMS program against a lateral trim solution computed at MSFC was 
made. The MSFC solution is for the case of zero net rolling and yawing moments, but, 
unlike in the TRIMS program the requirement of a zero net side force (i .e., Y = 0) is not 
imposed. Also the MSFC solution does not consider the deflection of the aileron. A special 
modification of the TRIMS program for including or disregarding the trim condition Y = 0 
and/or the aileron deflection was made and is described in Appendix C. Although the actual 
modification of the TRIMS program to eliminate the constraint Y = 0 is minor, it is based on 
a novel procedure derived in Appendix D. A comparison of the trim solutions computed by each 
program for (supposedly) the same trim problem showed that the deflection angles have about 
the same magnitude but are not equal . A more detailed discussion of the comporison including 
plots of the trim solutions is given in Appendix E, 
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FtguMi 4.6 txirafqt Trtwi of Spac* ShuttU 


CASE 1 


SVSTtM dynamics PAHAM£TE:RS 

— — — — — YAW BIAS TOHUUt = 3020000.0 

MOCL BIAS TOBUUt b 0,0 


PEfiFORMANCE CRITERION HAMAMLTEHS 


, — 

wn 

=3000. 00 

H2I 

a 

0.0 


H12 

=3000.00 

H22 

a 

0.0 

WEIGHTING 

Hl3 

=3000,00 

H23 

a 

0.0 

FACTORS 

WI4 

=3000.00 

H2A 

a 

0.0 


HIS 

=3000.00 • 

H2S 

a 

0.0 


HlG 

BAOOO.OO 

H26 

m 

0.0 


W17 

=3000.00 

H27 

m 

0.0 


TRIM deflection ANGLES 


TKAJo 

PT. 

FLIGHT 

time 

(1) 

(2) 

<3> 

delta 

(4) 

(5) 

(6) 

(7) 


’ 

V 

I 

25.0 

-11,50 

A. 08 

-10.35 

0.0 

0.0 

-7.55 

13.03 

• 



2 

40,0 

14.77 

4.99 

16.13 

0.0 

0.0 

4.28 

-13.35 


3 

50.0 

46,00 

3.32 

61,74 

0.0 

0.0 

31.07 

-24.96 * 


4 

OO.O 

100.03 

-4.25 

146.84 

0.0 

0.0 

81.55 

-39.35 


5 

65.0 

128.30 

-13.30 

229.11 

0.0 

0.0 

121.16 

-43.78 


6 

70.0 

95.05 

22, S8 

141.57 

u.o 

0.0 

18.55 

-68.60 • 


7 

75.0 

120.30 

8,95 

239.63 

0.0 

0.0 

20.78 

-70,80 


B 

BO.O 

136.46 

-2,99 

314.14 

0.0 

0.0 

35.18 

-80.27 


0 

00.0 

11.56 

22,01 

201.30 

0.0 

0.0 

52.10 

-67,40 . 


10 

100.0 

5.74 

-5.80 

-37.30 

0.0 

0.0 

-5.75 

27.50 • 


11 

110.0 

. 12.93 

-9,93 

-41,76 

0.0 

0.0 

6.74 

60.32 , 


12 

140.0 

UB.73 

-57.42 

-119.78 

0.0 

0.0 

58,39 

110.65 * 




TOP 

YAW 

PITCH 

YAH 

PITCH 

aileron 

RUDDER 





URUITEH 


SRM 

— 

> 



CASE 2 


SYSTEM DYNAMICS PARAMETERS 

tA« BIAS TORQUE » 2500000.0 

ROLL BIAS TORQUE = 700000.0 


performance criterion parameters 


— 

Wll 

=3000.00 

W21 « 

0.0 


W12 

=3000.00 

W22 « 

0.0 

WEIGHTING 

H13 

=3000.00 

W23 « 

0.0 

FACTORS 

W14 

=3000. go 

W24 » 

0.0 


Wl5 

sSOOO.UO 

W25 > 

0.0 


Wl6 

=4000.00 

W26 « 

0.0 


wl7 

=3000.00 

W27 • 

0.0 


TRIM deflection ANGLES 

— thaj. flight delta 



PT, 

time 

ID 

(21 

<31 

( 4 ) 

(5) 

(6> 

171 



1 

25.0 

-11,40 

5.19 

-19.19 

0.0 

0.0 

-5.80 

8.18 



2 

40.0 

14.20 

5.93 

15,00 

0.0 

0.0 

. 4.79 

-15,16 



3 

50,0 

47,25 

4.11 

61.77 

0.0 

0.0 

32.15 

-26.64 



4 

60.0 

100.54 

-3.70 

147.57 

0.0 

0.0 

83.70 

-40.56 



5 

65.0 

128.08 

-13,09 

230.70 

0,0 

0.0 

123.92 

-44,77 



6 

70.0 

94,54 

23,77 

141.11 

0.0 

0.0 

10.24 

-60.80 



7 

75.0 

128.06 

0,73 

240.02 

0.0 

0.0 

21.24 

-71.94 



8 

BO.O 

136,29 

-2.36 

315,70 

0.0 

0.0 

35.80 

-81.62 



0 

00.0 

9,60 

24.38 

205,51 

0.0 

0.0 

54.00 

-70.14 



10 

100.0 

U6B 

-3,15 

-34,26 

o.v 

0.0 

-4,23 

21.97 



11 

110,0 

7. 48 

-6.47 

-37,58 

0.0 

0.0 

6.20 

47,90 



12 

140.0 

91.36 

-44,12 

-09.66 

0,0 

0.0 

46,67 

87,91 


• « 



TOP 

VAil 

PITCH 

YAH 

PITCH 

AILCRON 

RUDDER 

• « P 


UNtflTER 
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Hflura 4.6 (ConHnu»cQ 


CASt i 

SYSTEK dynamics ♦*ANAMfcT£RS 


VAU BIAS TOKOOb s 0.0 

HULL BUS TUHUUE > BTOOOd.O 


performance criterion PAHAMbTEHS 






Wll 

sJOOO.OO 

W21 >. 

0.0 


W12 

=3000. 00 

W22 > 

0.0 

WEIGHTING 

WI3 

>3000,00 

W23 « 

0.0 

factors 

Wl4 

>3000,00 

W24 * 

0.0 


WI5 

>3000,00 

W25 • 

0.0 


Wlb 

>4000.00 

W26 * 

0.0 


W17 

>3000.00 

W27 > 

0.0 


trim deflection angles 


traj. flight 
Hr. time 


I 

£ 

3 

A 

b 

6 

r 

8 

10 

11 

12 


<l) 


<2) 


(31 


delta 


(SI 


(61 


(71 


as.o 

40.0 
so.o 

faO.O 

6S.U 

70.0 
7b, 0 

80.0 
BO.O 

100.0 

110.0 

140. 0 


2.46 
2b. 63 
60.27 
114.42 
142,74 

103.77 
139.00 

146.77 
12.67 
'^3.02 
-4,01 
-6.61 


2.S9 
2.99 
0.63 
-8.03 
-17.98 
21.98 
7.19 
-S.SO 
24.91 
2. OS 
2.16 
4.25 


3.16 
31.20 
81. 3b 
169.9b 
256.67 
156,4b 
259.74 
340,59 
240.96 
-4.11 
-3.75 
-0,63 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0,0 

0.0 

U.O 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


4.32 

6.6B 

39.90 

93.27 

134.91 

20.14 

22.50 

38.18 

62.35 

0.66 

0.88 

-2.07 


-10,51 

-22.68 

-32.28 

-45,32 

-48.87 

-75.25 

-76.72 

-87,15 

-80.79 

- 1.00 

-2,53 

-4.59 


TOP YAW 

< OrtGITER 


pitch yaw pitch aileron RUOOER 

>< 5RH 


CASE 4 


system dynamics parameters 

— — — — — YAtt BIAS TORQUE « -2500000.0 

moll bias Torque 700000.0 


performance CRITERION PARAMETERS 


— ... 

Wll 

>3000.00 

W21 > 

0.0 


WI2 

>3000.00 

«22 > 

0.0 

WEIGHTING 

Wl3 

>3000.00 

W23 « 

0.0 

factors 

Wl4 

>3000.00 

w24 > 

0.0 


wI5 

>3000.00 

«25- « 

0.0 


Wlb 

>4000.00 
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Figure 4.6 (Conlmoed^ 
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4.3 OPTIMUM FEEDBACK CONTROL AND PERFORMANCE 


The vector-matrix equations defining the linear stochastic optimum control problem 
and its solution are given in Section 3.4. These equations entail computation of the matrices 
F , M , K , P from the matrices A, B,C,E,G,H,V,W,Q,R defining the optimum 
control problem. In order to simplify the feedback design, the matrices F , M , K , P are 
partitioned as follows: 



A block diagram of the complete closed loop system in terms of the matrices listed above is 
given in Figure 4 7 . The lower half of the block diagram depicts the optimum feedback 
control system. 


For the lateral control of the Space Shuttle the state vector x , control vector u , 
bias vector z , and observation vector y are defined to be 
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The vehicle dynamics used in the optimum control design are for the first trajectory point 
(flight time - 25 sec). The elements of the matrices A , B , C were computed from the data 
In Appendix B using the equations in Section 4.1. It should be noted that the results in this 
section are for the case of no dorsal fins and the results in Section 4.2 include the effect of 
the dorsal fins. For the first trajectory point the difference in the two cases is minor. Since 
y is assumed to be a subvector of x , the elements of the observation matrix H are either 
0 or 1 where a value ofl in column i indicates that x^ is one of the measured quantities. 
The relative weighting matrices Q and R in the performance criterion (3.82) are selected by 
the control designer with the goal of optimizing the closed loop performance. The porticular 
approach adopted for this problem is to select Q and R of the form 


R = Diag [ u 


-2 

-2 -2 

^max ' 

<p / ^ / V 

^max ^max 

-2 

-2 

*^1 max 

' * * ' ' ^7 max 


-2 

max 


-2 -2 
^ max max 


The parameter a in the performance criterion (3.82) is varied until an acceptable "trade 
off" is achieved between the closed loop performance and the level of control effort. For 
this example ct = 1 and the maximum values of the state variables and conlirol deflections 
used in the performance criterion are listed below. 


max 




max 




max 


max 


max 


max 


- 10m 

= 0.01 rad 

= 0.01 rad 
= 5 m/sec 

= 0.1 rad/sec 

= 0.1 rad/sec 


u . = 0.2 rad 

jmax 


J=1 . ■ 


7 


The matrix V defining the state excitation noise spectral density is assumed to be diagonal 
with 
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V = 0 


i ^ 1, 2, 3 


'^44 ^^''max^ ^ 

''55=<°-'Pmax)' 

For this example the matrices G and E defining the noise model associated vvith the bias 
inputs are 

G = Diag CO.t , 0.1 X lo’^ , 0.1 X lo’'*] 

E = Diag 0.01 , - 0.01 , 0.01 ] 

The tth diagonal element of G is roughly equal to the max irnum value of squared. The 
negative diagonal elements in E provide a small amount of damping in the noise model which 
is required in order that the covariance matrix corresponding to the bias vector z does 
r>ot become infinite. The observation noise spectral density matrix is assumed to have the form 

w = Diag Cffy , er^ , er^ , ffp , af 1 

The standard deviation defines the level of noise associated with the measurement of y 
and the other standard deviations are similarly defined. By varying the standard deviations of 
the sensor noise as part of the design procedure,different Kalman filter designs are obtained. 
For the final Kalman filter design in this example 

a = 0. 1 y 
y ' max 

cr = 0. 1 ip 
<p max 

(T I = 0* 1 i/) 

ij> ^max 

O' = 0.01 p 

p max 

<j = 0.01 r 
r max 

The numerical values of these matrices defining the optimum control problem are given in 
Figure 4, 8 (a). 
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A computer program entitled Linear Systems Resign (LSD) was used to design the 

optimum feedback system. The LSD program solves the equations for the optimum control 

solution given in Section 3.4.1. The resulting matrices used to design 

the deterministic quadratic optimum control are shown in Figure 4.8(b). Similarly, the 

resulting matrices P, , / Pr, / l< / K used to design the Kalman filter are shown in 

I Z J X z 

Figure 4.8 (c). 

The performance achieved by the feedback control system was simulated for the 
different designs. The control deflections as a function of time were plotted and are shown 
Figures 4.9 and 4.10 . Note the SRM deflections 6 ^ and 6 ^ are not plotted since in the 
current investigatio'' it is assumed that the SRM are not gimballed. The dynamic response in 
Figures 4. 9 and 4.10 Is for the case where the vehicle starts from the trim condition for 

^ = 1 . 20 ° 

T = 3.02 X 10*^ N-m 
■■b 



to which correspond the deflection angles 
6 ^ = -6.63° 

62 = 7.35° 

63 = - 11-27° 

6 ,= 2.65“ 

6 - = - 12 . 11 ° 

The initial trim solution shown above is indicated by the straight lines in Figures 4.9 
and 4.10, The effect of a 2° step change in the sideslip angle causing an increase from 
^ = 1.20° to = 3.20° was simulated. The transient response curves show the performance 
of the control system in achieving the new trim solution. The curves in Figure 4.9 (a) are 
for the case of complete state feedback which assumes that the state of the process can be 
estimated perfectly (i.e. ,H=Z,V=W=0 ,x =x). This is not realistic but provides an 
upper bound on the performance as the estimation capability of the Kalman filter improves. 
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Observe in this case that the control deflections change discontinuaously due to a step chonge 
in jS . This does not occur when the Kalman filter is included. The curves for the remaining 
cases show the performance when different Kalman filter designs are used. The different 
designs correspond to different values of the W matrix as shown below 

V*/ rN. p 2 2 2 2 2 T 

W = Diag Cay , OTp , crj 


Figure 
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cr 
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<T 
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4.10b 
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•01 
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.1 

4.10c 
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.01 

.01 

.01 

.01 

4.11 
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.001 

.001 

.001 

.001 


The value of W in Figure 4.8 (a) and the matrices in Figure 4.8 (c) correspond to the Kalman 
filter design used In Figure 4.9.* 

In Section 3.4.1 a convergence property relating the optimum control approach and the 
trim control approach is given by (3.117). A demonstration of this property for the lateral 
control problem of the Space Shuttle is given below. Trajectory point number 1 occurring at 
25 seconds after launch is shown in which the roll and yaw bias torques due to misalignment 
of the sol id rocket motors are assumed to be 

roll bias torque = 3,02 x 10 (N-m) 

yaw bias torque = 0 

The control vectors , u^( “ ) , and u( «> ) obtained by the optimum control approoch 

have been computed in this case for three different values of a (1 0.1, 0.01) and are listed in 
Table 4.4. The computations were performed according to (3.94), (3.95), (3.98)-(3. 102) where 
the control weighting matrix R was chosen to be and where I denotes the identity matrix, 

R = 25X 

The trim solution or limiting solution for a - 0 was computed using the TRIMS program and is 
also listed in Table 4.4. An examination of Table 4.4 illustrates that the steady state control 
level u( » ) for the optimum control solution approaches the trim solution as cr approaches zero. 
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Figure 4.8 Numerical Volug of Matrices Used in Optimum Con^rol Design 


(g) Definitiori of Lineor Stochostic Optiinunr Con^fol Problem 
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Figura 4.8 Nmnericol VqIu« of Motric»i U»d in Optimum Control Da^ign, (CooHnuad) 


(b) DetefminUtic Quodrotic Optimum Control Daiign (Ccrwplata State Feedback) 
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Figure 4.9 Dynamic Response for a 2° Step in Sideslip Angle 
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Figure 4. 10, continued - 2 - 
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Table 4.6 Convergence of Optimum Steady State Control Level to 
Trim Control as Conlrol Weighting Decreases 
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a = OJ 


a = 0.01 


(trim solution) 
cr = 0 
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5. CONCLUSIONS AND FUTURE WORK 

Solutions to the trim problem can be efficiently calculated by the TRIMS program 
using the numerical methods described In Section^ 3. 1 . The results of this inyestigdtioh 
indicate that numerical solution by the Newton-Raphson method is preferable to the 
steepest descent method because it yields faster convergence and does not require the 
user to specify an iteration step size <r . If the initial guess of the solution used to start 
the Newton -Raphson method Is not in the region of convergence then the method may not 
converge or may converge to the wrong solution* In this case the steepest descent method 
should be used for the first few iterations to generate a good starting solution to the 
Newton-Raphson method. This hybrid method could be implemented in the TRIMS program 
with minor modifications. However, it appears that for most practical trim problems an 
initial guess of 6 = 0 is always in the region of convergence. 

For the linear trim problem, a diagonal weighting matrix .R in the quadratic per- 
formance criterion is sufficient in finding the "best" trim solution with respect to the 
limits on the deflection angles. Introducing nonzero values for the off-diagonal, elements 
of R complicates the selection of the performance criterion and does not , lead to a better 
trim solution than could be obtained by use of g diagonal matrix. Starting from the trim 
solution for a given diagonal R matrix, consider the problem of searching for a more 
desirable trim solution.. The penalty function method for varying the diagonal elements 
of R is a viable approach for improving the trim solution that Is easy to use. The 
penalty function method would be considerably facilitated if the computer computation 
of the trim solution Is performance in a conversational mode of operation rather than a 
batch mode. In the former case, the user can examine the trim solution and then immediately 
try a new R matrix. The process can be repeated in a single sitting as many times as is 
necessary. 

The lateral trim solution in Figure 4.6 indicates that the Space Shuttle configuration 
does not have sufficient control authority when the SRM engines are not gimballed, (The 
improvement in the trim solution obtained bygimballing the SRM engines is an area for future 
study which can be performed by the TRIMS program with minor modifications to the block 
data subroutine.) If the constraint of zero net side force (i.e., Y = 0) is eliminated and 
the vehicle is only trimmed in roll and yaw, the maximum control deflections decrease by 
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roughly an order-of~magnitude , In this case the trim solution is within the deflection 
limits. Hence, the control requirements increase significantly with the addition of the trim 
requirement Y = 0 . Maintaining Y = 0 is not as critical as zero net roll and yaw torques 
because angular errors are multiplied by the vehicle velocity in computing the displacement 
from the nominal trajectory. This suggests removing the trim condition Y = 0 entirely or 
replacing it by | Y | <€ . The value of ^ depends on how much side displacement error 
ia acceptable. By varying the weighting matrix R in the performance criterion with flight 
time rather than holding it constant, significant improvement in the trim solution might be 
achieved. The problem of realizing a trim solution for a time-varying R matrix must also 
be considered. 

Only the steady-state performance of the control system for bias inputs is considered 
in the trim calculation. The dynamic or transient response of the controls for fluctuating 
inputs must also be considered in the overall system design. For the nonlinear trim problem, 
the index of controllability defined in Section 3.3.2 is a quantitative measure for selecting 
the trim solution that results in the most controllable system with respect to the dynamic 
response problem. An integrol E proportional to the control energy is defined in Section 
3.3. 1 and is computed using the controllability Grammian W . Another measure for 
selecting the trim solution is given by the value of E . If the trim problem is linear, then 
the value of the controllability index or E does not vary with the trim solution. 

The basic question in studying the dynamic response problem is; ’What is the maximum 
deflection of each control for the possible fluctations in the disturbance inputs?" One 
possible approach to the dynamic response problem is to examine the values of E , where 
E denotes the energy expended by the 1th control to return the vehicle to trim. The values 

Z/ 

E. can be readily computed from the controllability Grammain W . Although this approach 

7 / 

has potential in gaining insight into dynamic response problem, it possesses two major 
limitations. First, there is no simple relationship between the energy E and the maximum 
value of the transient response curve showing the variation in the control deflection angle 
with time. Second, the control signal corresponding to E cannot be realized by a linear 

v 

feedback control system. The trim problem concerns only the static performance and can be 
studied without considering the detailed design of the feedback control system. The dynamic 
response problem, however, concerns the closed-loop transient response and is strongly 
dependent on the design of the feedback control system. 
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The most realistic method and possibly the only practical method for studying the 
dynamic response problem is to design the control system and simulate the closed-loop 
performance. The application of optimum control theory provides a method for the design of a 
linear feedback control system that can solve both the trim problem and the dynamic response 
problem. The correlation between the trim solution and the optimum control solution derived 
in Section 3.4.1 indicates how the solution to the trim problem can be used to select the 
proper control weighting R in performance criterion of the optimum control approach. 

This saves design time since the trim problem is easier to solve. A computer program 
entitled Linear System Design (LSD) was developed at Singer-Kearfott that is capable of 
computing the optimum feedback system and simulating the closed-loop performance. Since 
LSD is a conversational program with an automated plotting capability, many different 
designs can be studied efficiently. An example illustrating the use of LSD to design an 
optimum feedback system for the lateral control of the Space Shuttle during ascent is 
described in Section 4.3. It is recommended a more extensive design effort be pursued 
using the optimum control approach. 
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APPENDIX A VECTOR NOTATION AND DIFFERENTIATION 

In this appendix the notation used for handling differentiation with respect to vector 
quantities is reviewed for reference purposes. This notation is useful in describing the solution 
to the trim control problem. 

Let X and y denote an n dimensional and an m dimensional (column) vector, 
respectively. Further, let a denote a scalar function of x and y and let f denote a 
vector function of x and y where the dimension of f is p . 


f = f(x,y) 01 = C<x,y) 


w • 


• • 



H 




• 

v = 

- 

* 


* 

x 



n 


* m 



_ 


Differentiation of a vector by a scalar results in a (column) vector defined by 


X = dx/dt = 


dx^/dt 


dx /dt 


On the other hand, differentiation of a scalar by a vector results, in a row vector defined 


by 


^«/Sx = [Ba/^x^, Sfv/dx^/ • - • t Bof/^x^] 

The second partial of the scalar ot with respect to x and y 

Cif '^y.'by — by( da/ Sx)^ 
is an n by m matrix whose element is defined by 

( d^a/dxSy)^^ = d^a/dx^dy^ 

Differentiation of the vector function f with respect to the vector x is a p by n 
matrix whose tjfth element Is defined by 


(df/dx)^= df^/dx^ 
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Consider the scalar formed by the Inner product of f and a constant vector X of 
dimension p . The second partial of this scalar with respect to x and y 


a^(X'f)/3xSy = x'Of/axay) 

is an n by m matrix whose tjfh element is given by 



The quantity of 3f/Bx dy is a tensor whose ktjfh element is defined by 

Of/ax3y)^^^= 
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APPENDIX B PARAMETERS OF SPACE SHUTTLE DYNAMICS 


The equations defining the lateral -direction dynamics of the Space Shuttle during 
ascent through the atmosphere were derived in Section 4,1 . The parameters required to 
compute the matrix coefficients in the linear equations of motion (4. 10) are given in this 
appendix. The list of parameters appearing below indicates the parameter symbol, value, 
units, and a brief description. The data is given for twelve different points or flight times 
along the ascent trajectory and was furnished by Dr. S. Winder of MSFC. 


In the column labeled VALUE, there appears either the numerical value or the word 
"table" or is left blank. The word "table" denotes that the numerical value varies with 
flight time and the twelve different values are listed in the tables at the end of this appendix. 
A blank denotes that the value of the parameter has not been specified. The unspecified 


parameters are the location and thrust of the SRM engines and the stability derivative C^ 
Most likely C is small and is assumed to be zero In this investigation. It is 
further assumed ^ that the SRM engines are not gimballed but the provision for including 



the SRM engine deflections is incorporated into the equations. 


The stability derivatives C ,C ,C ,C ,C were not included in the data 
' tp np yr ^ nr 

furnished by MSFC. Their values listed below are rough estimates based on the vehicle 
configuration. These stability derivatives are not used in computing the trim solution but 
are required for the study of dynamic response. 
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LATERAL TRIM PARAMETERS 

SYMBOL VALUE 

UNITS DESCRIPTION 

X, 0 

m ) 


^1 

0 

m ' 

/ 

^ X, y, 2 positions of (top, orbiter) engine 1 


0 

m 

m ' 

) 

CN 

X 

0 

m 1 

} 

^2 

1,346 

m ' 

^ X, y, z positions of (right orbiter) engine 2 

^2 

- 6.68 

^ ( 

• > 


’'3 

0 

m 

1 

) 

^3 

- 1.346 

m ' 

^ Xf y, z positions of (left orbiter) engine 3 

^3 

- 6.68 



>'4 


m 


^4 


m 

1 

^ X, y, z, positions of (right- SRM) engine 4 

^4 


m 1 

1 

>‘5 


m ^ 

) ' 

^5 


m 

^ X, y, z position of (left SRM) engine 5 

^5 


m 

^ ■ 

X 

eg 

table 

m j 


^cg 

0 

rn 

^ X, y, z position of center of gravity 

z 

eg 

table 

m 


X 

mrp 

21.6 

m 

/ 

y 

'mrp 

0 

m 

\ X, y, z position of moment reference point 

z 

mrp 

- 1.47 

m 
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q 

table 

/ 2 

new ./m 

dynamic pressure 

s 

317.73 

2 

m 

reference area 

b 

28.322 

m 

reference length 

V 

table 

m/sec 

velocity of the vehicle relative to the air 

V 

fable 

m/sec 

side component of V (side wind velocity) 

y 

F 

table 

New. 

thrust per orbiter engine 

’'SRM 


New. 

thrust per SRM engine 


table 

- 

stability derivative 


table 

- 

stability derivative 


table 

- 

stability derivative 


table 

- 

change in due to dorsal fins 


table 

- 

change in C _ due to dorsal fins 

^^'"n^^AFT 

. table 

- 

change in dorsal fin 

^^‘"n^FORWARD 

table 

-■ 

change In C _ due to forward dorsal fin 
nd 

C 



stability derivative 

. ^ea 

c. 

table 


stability derivative 

fia 

C 

table 


stability derivative 

fia 

C 

table 

- 

stability derivative 

’'■5T 

c 

table 


stability derivative 

fir 

C 

table 

— 

stability derivative 


fir 


'IP 


np 




--01 

- 0.03 


stability derivative 
stability derivative 
stability derivative 


0 . 
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0,022 

- 

stability derivative 

c 

rvr 

-o,n 

- 

stability derivative 

c 

20. 

m 

length of mean aerodynamic cord 

m 

table 

Kg 

vehicle mass 

i 

I 

X 

table 

Kg-m^ 

vehicle moment of inertia about x axis 

I 

y 

table 

Kg-m^ 

vehicle moment of inertia about y axis 

j 

z 

table 

Kg-m^ 

vehicle moment of inertia about z axis 

9 

table 

/ 2 
m/sec 

acceleration of gravity 

cos 0 
o 

table 

- 

cosine of nominal pitch angle 

sm 0 
o 

table 

- 

sine of nominal pitch angle 

Q 

o 

table 

rad/sec 

nominal pitch rate 

U 

0 

table 

m/sec 

nominal velocity along x axis 

W 

o 

table 

m/sec 

nominal velocity along z axis 

max 

30 

deg 

maximum allowable rocket engine deflection (1=1/ 

max 

table 

deg 

maximum allowable aileron deflection 

^ max 

table 

deg 

maximum allowable rudder deflection 


0 

2 

m 

reference area for drag induced by aileron control 

67 

0 

2 

reference area for drag induced by rudder control 

^D6 

0 


drag coefficient for aileron control 

^D 7 

0 

- 

drag coefficient for rudder control 


- 0.1 

- 

stability derivative 

C 

np 

- 0.03 

- 

stability derivative 

C 

yr 

0.0 

- 

stability derivative 

S 

0.022 

- 

stability derivative 

C 

-0.11 


stability derivative 



DATA: 



Stability Derivatives 
(all data/radian) 

\ % 


.273 

- .510 

- 1.66 

- .283 

.302 

.265 

- .489 

- 1.68 

- .285 

.315 

.259 

- .473 

- 1.70 

- .286 

.325 

.215 

- ,388 

- 1.83 

- .291 

.404 

.181 

- ,310 i 

- 1.99 

- .298 

.468 

.173 

- .345 

- 2.05 

- .326 

.460 

.206 

- .340 

- 1.97 

- .384 

.344 

.186 

- .254 

- 1.92 

- .356 

.266 

J 05 

- .137 

- 1.93 

- .299 

.238 

.055 

- .105 

- 2.03 

- .246 

.269 

.0406 

- .077 

- 1.98 

- .196 

.207 

.0286 

- .061 

- 1.60 

- . 122 

. 

- .0284 


DATA: Vehicle Parameters 


flight 

time 

(sec) 

m 

(Kg) 

I 

^ 2 
(Kg-m ) 

I 

y 2 

(KgV) 

I 

"" 2 
(Kg-m ) 

X 

eg 

(m) 

z 

eg 

(m). 

F 

(New.) 

25 

.218E+7 

.953E+8 

.526E-t9 

1 

.591 E+9 

23.345 

1 

- 1.58 

1.650E46 

40 

.201 E+7 

.856E+8 

.490E+9 

.547E49 

23.42 

- 1.5847 

1.760E46 

50 

.190E+7 

.794E46 

.468E49 

.519E49 

23.47 

- 1.5914 

1.825E46 

60 

.179E+7 

.733E48 

.445E+9 

.491 E49 

23.52 

- 1 .5953 

1.885E46 

65 

.174E+7 

.702E4B 

.434E+9 

.478E49 

23.545 

- 1 .5979 

1.920E46 

70 

. 169E+7 

.671 E+8 

.423E+9^ 

.464E+9 

23.57 

- 1.60 

1.940E^6 

75 

.160E+7 

.629E46 

.383E49 

.420E+9 

24.13 

- 1 .4626 

. 1.970E-WS 

80 

. 154E+7 

.606E+8 

.372E+9 

.405E49 

24.18 

-1.455 

1.980E46 

90 

.144E+7 

.559E48 

.348E+9 

.375E49 

24.33 

- 1.440 

2.025E46 

100 

.133E+7 

.512E+8 

.326E49 

.346E+9 

24.535 

- 1 .4327 

2.040E46 

no 1 

.122E+7 

.466E-t8 

.303E-+9 1 

.317E49 

24.74 

- 1.4255 

2.060E46 

140 

.914E46 

.329E-+8 

.234E+9 

1 

.228E49 

25.62 

- 1.400 


2.070E46 





DATA: Trajectory Parameters 


flight 

time 

(sec) 

V 

(m/sec) 

V 

y 

(m/sec) 

g 

(m/sec^) 

cos 0 
0 

sin 0 
o 

Q 

o 

(rad/sec) 

U 

o 

(m/sec) 

W 

o 

( m/sec) 

q 

(New/m^) 

25 

95,4 

2.0 

9.8 

- .012 

1.0 

- .203E-1 

95.4 

.279E40 

.482E+4 

40 

150, 

9.0 

9.8 

.050 

.999 

- .512E-2 

149. 

. 149E+3 

.987E44 

1 

50 

190. 

15.0 

9,79 

.125 

.992 

.450E-2 

186. 

.186E+3 

.134E+5 

60 

241. 

24.0 

9.78 

.222 

: .975 

.112E-2 

232. 

.232E+3 

.174E+5 

65 

272. 

29.0 

9.78 

.274 

.962 

- .158E-1 

257. 

.257E+3 

.194E45 

70 

305. 

34.0 

9.78 

.329 

.944 

.687E-1 

283. 

.283E+3 

.212E+5 

75 

343. 

40.0 

9.77 

.384 

.923 

- . 103E-+0 

310. 

.310E+3 

.226E+5 

80 

385. 

44.0 

9,77 

.449 

.893 

.412E-1 

337. 

.337E+3 

.233E+5 

90 

486. 

30.0 

9.76 

.566 ' 

.824 

- .227E-2 

392. 

.392E+3 

.217E+5 

TOO 

612. 

0.0 

9.74 

,664 

.748 

.301E-3 

445, 

.445E+3 

.165E+5 

no 

768. 

0.0 

9,73 

.681 

.732 

- .851E-2 

498. 

.498E+3 

.117E+5 

140 

1520. 

0.0 

9.68 

.874 

.486 

- .lOOE-2 

673. 

.673E+3 

.231E44 


DATA: Deflection Limits for Aerodynamic Surface Controls and Change in Stability Derivatives Due to Dorsal Fins 


flight 

time 

(sec) 


7 max 

rudder hinge 
moment limit 
(deg) 


6 max 

aileron hinge 
moment limit 
(deg) 



^all data / degrees) 

^^^niS^AFT FORWARD 


25 

nohinge limit 

1 nohinge limit 

- .011 

.0031 

.0064 

o 

o 

1 

40 

42.0 

! 71.8 

- .012 

.0032 

.0067 


50 

30.8 

52,6 

- .013 

.0033 

,0074 


60 

23.5 

40.0 

- .015 

.0036 

.0085 


65 

14.7 

25.1 

- ,016 

.0038 

.0094 

- .006 

70 

8.19 

14.1 

r .017 

.0042 

.0104 

- .006 

75 

5.54 

9.47 

-.0165 

.0042 

.01 

-.0058 

80 

5,23 

8.91 

- .014 

.0035 

.0088 

- .005 

90 

6.27 

10.69 

- .0105 

.0027 

.0075 

- .0044 

100 

10.23 

17.5 

- .008 

.0017 

.005 

- .0028 

no 

19.67 

33.64 

- .006 

.0014 

.004 

-.0022 

140 

no hinge limit 

no hinge limit 

- .004 

.001 

.0028 

- .0015 

* hard limits 

±30 

j 

40 up- 15 down 

n 


* Hard deflection angle limit is used when less than 
hinge moment limit. 




APPENDIX C TRIMS COMPUTER PROGRAM 


1. PROGRAM USAGE 
Input 

The input data to the TRIMS program consists of punched cards. The data deck 
is divided into cases where for example each case computes the trim solution for 
different values of roll bias torque. There are seven punched cards per case with the 
first card containing the case title and the last card indicating whether another 
case follows or whether this is the last case to be run. A description of the information 
and format for punching these seven data cards per case is given in Table 1. A sample 
of an input data deck for a single case run is shown in Figure 1. 

Output 

The computer printout from the TRIMS program is a single page per case. The 
printout resulting from the data deck in Figure 1 is shown in Figure 2. The first part 
of the printout lists the information contained on the data cards and used to compute 
the trim solution. The trim solution is printed in a convenient tabular form with each 
row listing the seven trim angles in degrees for a particular flight time. The number 
of iterations required to compute the trim solution at each trafectory point is also 
indicated. 

Options 

Special options have been added to the program since the original development 
date of February, 1973, The purpose of these options is described in Table 2 including 
the modifications to the input data required to exercise these options. 


U2 



TABU 1; TRIMS PROGRAM INPUT DATA 
CARO COLUMNS VARIABLE FORMAT DESCRIPTION 


TITLE CARD 


1 1-72 LINE 72A1 

CONTROL CARD 


1 

1-5 

IGRAD 

I 5 

1 

11-20 

EPS 

El 0.3 

1 

21-30 

STEP 

El 0.3 

TRAJECTORY CARD 



1 

1-60 

JPT 

12j5 

CARD CONTAINING BIAS 

: TORQUES 


1 

1-10 

YBT 

El 0.0 

1 

11-20 

RBT 

ElO.O 


CARDS CONTAINING WEIGHTING FACTORS 


Descriptive case title, ^ 

= 1, use 1st order gradient method; 

= 2, use 2nd order gradient method. 

Upper bound used in the convergence criterion. 

Step size used in the 1st order grodient method; 
leave blank if 2nd order gradient method is used. 

If trajectory point no. k, lc=l, 12, is to be 
used then punch a 1 in column 5k ; other- 
wise punch a 0 in column 5k. 

Yaw bias torque. 

Roll bias torque. 


1 

1-70 

W1 

7E10.0 

Seven weighting factors in performance criterion 
for adjusting maximum deflection angles. 

2 

1-70 

W2 

7E10.0 

Seven weighting factor in performance criterion 
for adjusting aerodynamic (drag) and thrust 
losses due to trim. 

CASE PARTITION CARD 




1 

1 

IG0 

T5 • 

r= 1, another case follows 
L= 2, last case. 


FIGURE 1; EXAMPLE OF INPUT TO TRIMS COMPUTER PROGRAM 


STUOY ftF 


TPIM 



NO 

SPM RIAS 


? 

0.0001 

0. 






1 

11 

1 1 

1 

1 1 

1 

1 1 

1 

ft. 

ft. 







3000. 

3000. 

3000 

• 

3000. 

3000. 

3000. 

3000 

0. 

ft. 

ft. 


(1. 

0. 

0. 

0. 


? 
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FIGURE 2: EXAMPLE OF OUTPUT FROM TRIMS COMPUTER PROGRAM 


CASF 1 STUDY OF LATERAL TRIM FOR SPACE SHUTTLE NO SRM BIAS 


COMPUTATION CONTROL PARAMFTEPS 


use 3ND order r^OAQIENT METHOD 


trajectory points 


UPPER BOUND USED IN THE CONVERBENCE CRITERION = O.lOOE-03 


1 


1 


SYSTEM dynamics PAPftMETERS 

YA^ BIAS TORQUE ^ 0.0 

ROLL BIAS TORQUE = 0.0 


perrormamce criterion parameters 



Wll 

=3000,00 

W21 

= 

0.0 


Wl? 

=3000.00 

W22 

= 

0.0 

wFISHTING 

W13 

=3000.00 

W23 

= 

0.0 

factors 

WlA 

=3000.00 

W24 

= 

0.0 


wis 

=3000,00 

W25 

= 

0.0 


Wlh 

=3000,00 

W26 

= • 

0.0 


W17 

=3000.00 

W27 

= 

0.0 


trim deflection anglfs 


TRA J, 
PT. 

FLIGHT 

time 

( 1 ) 

<2 > 

(3) 

delta 

(^) 

(5) 

(61 

(7) 

NO, OF 

iterations 

1 

26.0 

0,09 

-0,00 

0,13 

0,0 

0.0 

0.10 

-0.10 

« 

. 1 

? 

40.0 

0.4S 

-0.04 

0.59 

0.0 

0.0 

0.14 

“0.26 

. 1 

3 

50.0 

0.87 

“0.13 

1.24 

0.0 

0.0 

0.61 

-0,26 

. 1 

4 

60,0 

1.5H 

-0,31 

2.40 

0.0 

0.0 

1.43 

“0.36 

. 1 

5 

65.0 

2,04 

-0,48 

3.50 

0.0 

0.0 

1.86 

“0.39 

. 1 

h 

70.0 

1.71 

0,02 

2.57 

0.0 

0.0 

0.25 

-o,e« 

, 1 

7 

75.0 

2.26 

-0,24 

4.03 

0.0 

0.0 

0.39 

“0.89 

. 1 

8 

80,0 

2.42 

-0.39 

5,15 

0.0 

0.0 

0,78 

-1.03 

. 1 

9 

90.0 

0,54 

0.12 

3,48 

0.0 

0.0 

1,31 

-0.94 

* 1 

10 

100.0 

“0,00 

-0,00 

0.00 

0,0 

0.0 

0.00 

-0.00 

. - 1 

1 1 

110.0 

0,00 

0,00 

-0,00 

0.0 

0.0 

0.00 

0.00 

. 0 

12 

140.0 

0,00 

0,00 

-0.00 

0.0 

0.0 

0.00 

0.00 

• 0 


TOP yaw pitch yaw pitch AILERON RUDDER 

< — — ORRITER SRM , > 



TABLE 2: PROGRAM OPTIONS 


Option 1 - The program has the capability of disregarding the first trim equality constraint. 
This equation corresponds to the trim condition of zero net force in the y-^irection. To 
exercise this option change the nonzero values of JPT on the trajectory data card from 
positive numbers to negative numbers. 

Option 2 - The program has the capability of computing the trim solution for the case 
where the aileron is not used. To exercise this option change the nonzero values of JPT 
on the trajectory data-card from a magnitude of 1 to a magnitude of 2 (i.e., replace 
1 by 2 and replace - 1 by - 2). 

Option 3 - The program has the capability of replacing the performance criterion stored 
internally in the program with the quadratic performance criterion 

r(6) = (a,/e,)^ + . . . +(Sy/cy)^ 

where c^ , . . . , c^ are seven constants specified by the user at execution time. To 
exercise this option replace the fourteen weighting factors in the input data with the 
values 


Wl(i) = -c^ 
t 

W2(t) = 0. 
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2. PROGRAM DESCRIPTION 


TRIMS IS a FORTRAN IV computer program composed of a single mam or 
executive routine and many subroutines. The program subroutines may be viewed 
as divided into two main groups. The first group is comprised of the main routine, 
entitled TRIMS, plus seven basic subroutines which form the heart of the program. 

These are listed in Table 3 together with a brief description of their function. The 
second group contains the utility subroutines which perform o specific matrix operation 
such as invert a matrix or print out a matrix. There are thirteen of these subroutines 
which are listed in Table 4. With the exception of GMSYMM, all of the utility 
subroutines are found In the IBM Scientific Subroutine Package . 

In addition to the calling lists, the transfer of information into and out from 
the subroutines is achieved by means of five named C0MM0NS. Their names are 
listed in Table 5 together with a brief functional description. The innerconnection 
between the main routine, the seven basic subroutines, and the five named C0MM0 NS 
summarizing where each is used is shown in Tabled . The variables in each of the 
named C0MM0NS are listed and defined in Table 7. The other variables in the 
program not in a named C0MM0N are listed in Table 8. 

In the following pages the FORTRAN source listing of each subroutine is given. 
The beginning of each listing contoins comment cards describing the subroutine which 
includes the purpose, input variables, output variables, and the subroutines called. 

Flow diagrams are also given for each of the subroutines with the exception of the 
IBM SSP subroutines. 


System/360 Scientific Subroutine Package, Version iZJ, Programmer’s Manual, 
I BM publication GH20-0205-4, Fifthe edition, August 1970. 
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TABLE 3 : MAIN ROUTINE AND BASIC SUBROUTINES 


TRIM 

BLOCK 

INPUT 

OUTPUT 

GRAD1 

GRAD2 

SYSTEM 


COST 


GYSYMM 

MCPY 

MSTR 

L0C 

GMSUB 

GMPRD 

GMTRA 

MPRD 

CCUT 

MINV 

SINV 

MFSD 

MX0UT 


main roufine controlling the basic computotional steps. 

block data subroutine for storing data internally in the program, 

subroutine used to read in and print out the input dato. 

subroutine used to print out the results of the program. 

subroutine for computing the deflection angles using the 1st order 
gradient method. 

subroutine for computing the deflection angles using the 2nd order 
gradient method, 

subroutine containing the equations defining the system dynamics 
and the corresponding equations for evaluating the derivatives 
required by the gradient methods, 

subroutine containing the equations defining the performance 
criterion and the corresponding derivatives. 


TABLE4 : UTILITY SUBROUTINES 

symmetrize a matrix 
matrix copy 

storage conversion of a matrix 
location in compressed -stored matrix 
subtract two general matrices 
product of two generol matrices 
transpose of a general matrix 
matrix product 
partition a matrix by column 
matrix inversion 

invert a symmetric positive definite matrix 

triangular factorization of a symmetric positive definite motrix 

print o matrix 
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/C0N/ 

/array/ 

/TRAJ/ 

/SYST/ 

/PERF/ 


TABLE 6: 


where 

used 

requires 

TRIMS 

BLOCK 

INPUT 

0UTPUT 

GRADl 

GRAD2 

SYSTEM 

C0ST 


TABLES : NAMED C0MM0NS 

dimension and accuracy parameters 

values of trim equation, performance criterion, and 
their derivatives 

trajectory information 

data derived from the space shuttle configuration for 
computing the system dynamics and trim equation 

data used to compute the performance criterion 


INNERCONNECTION OF SUBROUTINES AND NAMED COMMONS 


1/1 


a. 


Z) 

Q- 





<N 


C0ST 





a 

S 

o 

S 

is 

1— 

t/1 

z 

•O 

oc 

— \ 
s 

1— 

to 

>- 

O 

o 

> 

t/1 


u 

X 

< 

1— 

to 




1 





X 

X 


1 

X 







1 

1 

X 


X 

X 




1 

1 

1 

X 


X 

X 



X 

1 

xl 

X 

X 

X 




X 

X I 

X 

X 

X 





1 

1 


X 


X 


X 


O' 

LU 

O- 


X 

X 


X 
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TABLE 7: Voriables In Named C0MM0N 


r 


Program 

Symbol 

Dimension 

Symbol 

/C0N/ 

M 

* • « 

m 

NS 

• ♦ * 

n 

KMAX 

* e « 

K 

max 

EPSO 


% 

MPT 

• • • 

« • ♦ 

/ARRAY/ 

AV 

6 

a 

av 

6 

|(6) 

BM 

60 


BT 

6,60 

d^b/^ 

RS 

« • » 

r 

RV 

10 

dr/ 

RM 

100 

6^/66^ 

/TRAJ/ 

JPT 

12 

• « e 

TF 

12 

• « * 

/SYST/ 

YBT 


0 0 0 

RBT 


0 0 0 

S 


S 

BREF 


^ref 

X1,Y1,Z1 



X2,Y2,Z2 



X3,Y3,Z3 



X4,Y4,Z4 



X5,Y5,Z5 



XMRP 


X 

mrp 

YMRP 


y 

mrp 

ZMRP 

• e • 

z 

mrp 

XCG 

12 

X 

eg 

ZCG 

12 

eg 

11 


Explanation 

Number of trim angles. 

Number of trim equations. 

Maxirhum number of iterations ol lowed. 

Relative tolerance used in subroutine SINV. 
Maximum number of trajectory points a I lowed. 

Constant terms in trim equations. 

Terms in trim equations varying with trim angles 

First derivatives of trim equations. 

Second derivative of trim equations. 
Performance criterion. 

First derivative of performance criterion. 

Second derivative of performance criterion. 

Index vector determining which trajectory 
points to use (see program input data). 

Flight times corresponding to the different 
possible trajectory points. 

Yaw bias torque (see program input data). 

Roll bios torque (see program input data). 
Reference area. 

Reference length. 

Coordinates of (top orbiter) engine 1. 
Coordinates of (right orbiter) engine 2. > 

Coordinates of (left orbiter) engine 3. 
Coordinates of (right SRM) engine 4. 
Coordinates of (left SRM) engine 5. 

' ■ 

» Coordinates of moment reference point. 

{ Coordinates of center of gravity 
eg 


r 



TABLE/: 

Variables 

In Named C0MM0N, Confinued 

Program 

Symbol 

Dimension 

Symbol 

Explanation 

/SYST/, Confinued 

Q 

12 

q 

Dynamic pressure. 

V 

12 

V 

Vehicle velocity relative to air. 

VY 

12 

V 

y 

Side wind velocity. 

F 

12 

F 

Thrust per orblter engine. 

FSRM 

12 

''SRM 

Thrust per SRM engine. 

CYB 

12 


Stability derivative. 

CLB 

12 

C 

Stability derivative. 

CNB 

12 

% 

Stability derivative. 

DCYB 

12 


Change In C due to dorsal fins. 
y/3 

Change In C due to dorsal fins. 

DCLB 

12 

DCNBA 

12 

(aC )aft 

FORWARD 

0 

Change In C due to aft dorsal fin. 

Change In C due to forward dorsal fin. 
r»j3 

DCNBF 

12 

CYA 

12 

C 

y6a 

Stability derivative. 

CLA 

12 

C 

Stability derivative. 

CNA 

12 

C 

nsa 

Stability derivative. 

CYR 

12 

C 

y6r 

Stability derivative. 

CLR 

12 

c 

n§r 

Stability derivative 

CNR 

12 

Stability derivative. 

/PERF/ 

W1 

7 

w, 

W 

Vector of relative weighting factors 

W2 

7 

(see program Input data). 

Vector of relative weighting factors 

DAMAX 

12 

Z 

6 

(see program input data). 

Maximum deflection angle allowed for aileron. 

DRMAX 

12 

a max 

6 

Maximum deflection angle allowed for rudder. 

QQ 

12 

rmax 

q 

Dynamic pressure. 

DMAX 

• • • 


Maximum deflection angle allowed for 

SA 

« * * 

S 

orblter rocket engines. 

Reference area corresponding to the drag 

SR 

* ♦ « 

Q 

s 

induced by the aileron. 

Reference area corresponding to the drag 



r 

induced by the rudder. 

CDA 

• • « 

^Da . 

Coefficient of drag corresponding to the aileron. 

CDR 

» • « 

^Dr 

Coefficient of drag corresponding to the rudder. 
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TABLE 8 * Variables not in Named C0MM0N 


Pipgrom 

Symbol 

JCASE 

7 GO 

I 

L 

7GRAD 

K 

EPS 

TIME 

STEP 

DELTA 

LAMDA 

J 

DET 

MNS 

N0RM 

RU 

RX 

X 

X 

BX 

DU 

BU 


J7 

M2 

7ER 

HL 

R 


Dimension Symbol 


# # # 


10 

6 


k 

€ 

t 

a 

d 

X 


10 r 

u 

10 r 

X 

10 X 

10 

60 B 

X 

10 

60 B 

u 


*04 004 

0 4 0 0 4' 0 


10 

’00 >’66 


Explanation 

Number of the current case. 

Index controlling sequence of cases. 

Do loop index. 

Number of the trajectory point. 

Order of the gradient method to be used. 

Number of the iteration. 

Convergence bound in gradient methods. 

Flight time of the current trajectory point. 

Iteration step size used in first order grodient method 
Vector of trim angles. 

Vector of Lagrange multipliers. 

Do loop index. 

Determinant of a matrix 

Difference between number of trim angles and trim 
equations 

Quantity for determining trim solution accuracy. 
Subvector of dr/ B 6. 

Subvector of B r/d 6. 

Subvector of 6 (subroutine GRAD 1). 

Dummy vector (subroutine GRAD2). 

Square nonsingular submatrix of Bb/d6. 

Correction to subvector u of 6 . 

Submatrix of B b/ d . 

Matrix element index 
= m (m + 1 )/2, 

1 

Index used to indicote errors in inverting 
a positive definite matrix. 

Derivative of hamiltonian with respect to X . 

Second derivative of hamiltonian with respect to 6 * 
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TABLE 8 : Variables nof in Nomed C0MM0N (Confinued) 



Dimension 

Symbol 

Explanation 

Y 

10 

. . . 

Dummy vector 


DEL 

10 

A6 

Correction to 6» 


BR 

60 


Matrix product. 


LAM 

6 

BR'^B' 

Correction to X . 


BRB 

36 

Motrix product. 


HD 

10 


Derivative of Hamiltonian with respect to 

D 

60 

. • . 

Dummy matrix. 


B 

60 

^6 

Mixed second derivative of Hamiltonian. 

CYBCG 

« * * 


Stability derivative 

C about eg . 

y 

CNRCG 

# • « 


Stobility derivative 

C about eg . 
n 

r 

CLBCG 



Stability derivative 

C about eg . 

CNBCG 


. . . 

Stability derivative 

C about eg . 

n 

$ 

CLACG 

# • • 

# • * 

Stability derivative 

C about eg . 

^0 

CNACG 

• • • 

0 0 0 

Stability derivative 

C about eg . 

n 

a 

CLRCG 

• « • 

0 0 0 

stability derivative 

C about eg . 

Cl 

• « • 

« * . 

Cos 18°. 


C2 

0 0 0 

0 0 0 

Cos 12°. 


C3 

*00 

• 00 

Cos 3.5° 


C4 

0 0 0 

0 0 0 

Cos 15°. 


SI 

0 0 0 

'000 

Sin 18°. 


S2 

0 0 0 

0 0 0 

Sin 12°, 


S3 

0 0 0 

0 0 0 

Sin 3.5°. 


S4 

0 0 0 

0 0 0 

Sin 15°. 
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TABLE 8 • Variables not in Named C0 MM0N (Continued) 


Program 

Symbol 

Dimension 

Symbol 

Explanation 

IJ 

• # • 

« « * 

Matrix element index. 

QS 

e e « 

• • • 

Product qS . 

QSB 

• * e 

• • • 

Product qSb^^^ . 

RAD 

e • * 

♦ • « 

Conversion factor from radians to degrees. 

BETA 

* e • 


Side slip angle. 

II 

« « « 

• * # 

Vector element index. 

C 

« • ♦ 

e • • 

Dummy vector. 
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C TRM 0010 

C TRM 0020 

C TRIMS COMPUTER PROGRAM *** TRM 0030 

C TRM 0040 

C DEVELOPED BY: M, HUTTON THE SINGER CO. FEBRUARY 1973 TRM 0050 

C TRM 0060 

C TRM 0070 

C TRM 0060 

C PURPOSE main ROUTINE FOR EXECUTION OF COMPUTIONS OF LATERAL TRM 0090 

c trim angles for space shuttle. TRM 0100 

C TRM 0110 

C inputs (SEE SURPOUTINE INPUT). TRM 0120 

C TRM 0130 

C TRM 0140 

C OUTPUTS (SEE SUBROUTINE OUTPUT). TRM 0150 

C TRM 0160 

C TRM 0170 

C SUBROUTINES CALLED INPUT t GRAOl » GPA02 ♦ OUTPUT . TRM 0180 

C TRM 0190 

C TRM 0200 

C TRM 0210 

C TRM 0220 

G TRM 0230 

C TRM 0240 

REAL LAMDA TRM 0250 

DIMENSION DELTA(IO) » LAM0A(6) TRM 0260 

COMMON /CON/ M * NS f KMAX * FPSO f MPT TRM 0270 

TRM 0280 
TRM 0290 

INITIALIZATION TRM 0300 

TCASE = 1 TRM 0310 

DO 10 I=1*M TRM 0320 

10 DELTAU) = 0, TRM 0330 

DO 20 T=1*NS TRM 0340 

20 LAMDA(I) = 0. TRM 0350 

TRM 0360 

<MM» enter INPUT DATA TRM 0370 

30 CALL . INPUT (IGRAO^EPS^STEP^ IGO« ICASE) TRM 0380 

TRM 0390 

COMPUTE TRIM SOLUTION FOR EACH OF THE SELECTED POINTS ALONG TRM 0400 

THE TRAJECTORY TRM 0410 

00 60 L=1 *MPT TRM 0420 

TRM 0430 

DETERMINE COMPUTATIONAL METHOD TO BE USED TRM 0440 

lE(IGRAD-l) 50*40?50 TRM 0450 

TRM 0460 

COMPUTE TRIM SOLUTION USING 1ST ORDER GRADIENT TRM 0470 

40 CALL GRADI (K*L*TlMEtDELTA*LAMOA*lGRAD«EPS*5TEP) TRM 0480 

GO TO 60 TRM 0490 

TRM 0500 

*** COMPUTE TRIM SOLUTION USING 2ND ORDER GRADIENT TRM 0510 

50 CALL GRA02(K*L»TlME*DELTA*LAMnA*IGRAD»EPS) TRM 0520 

TRM 0530 

PRINT RESULTS TRM 0540 

60 CALL OUTPUT (K«L*M.TIME*OELTa*MPT) TRM 0550 

TRM 0560 

*** TEST IF END OE COMPUTER RUN TRM 0570 

IFdGO-l) 80*70*80 TRM 0580 

C TRM 0590 
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c r,0 TO THP NrXT CASF 

70 ICASF = TCASE ♦ 1 
on TO 30 
C 

«0 CALL EXIT 
END 


TRM 

TRM 

TPM 

TPM 

TPM 

TPM 


0600 

0610 

06?0 

0630 

0640 

0650 



JCASE = 
JCASE + 1 










O <J O O U O U O U U U O u u o oo 


I* 


c 

c 

c 

C 


SURPROGR&M BLOCK 


PURPOSE 


BLOCK data subroutine FOR STORING DATA INTERNALLY 
PROGRAM. THIS IS THE AERODYNAMIC DATA REQUIRED TO 
COEFFICIENTS OF LATERAL TRIM EQUATIONS. 


COMMON BLOCKS 


/CON/ 


/TRAJ/ • /SYST/ * /PERF/ 


BLOCK DATA 


COMMON /CON/ M 
COMMON /TRAJ/ 
COMMON /SYST/ 


« NS « KMAX « FPSO 
JPT(12) . TF(12) 


MPT 


XCG(12) ♦ 
FU2) * 
OCYB< 12) * 
CLAn2)* 
COMMON /PERF/ 


YBT* RBT* 

S* BPFFf 

ZCGU2) 9 
FSRM(12) ♦ 
DCLP(12> • 
CNA < 12) f 
(7) * W2(7) * 
DMAXf SA* SR» 


XI* 

Yl* 

Zl* 


X3* 

Y3* 

Z3* 


X2* 

Y2* 

Z2* 
Q(12) * 
CYB<12) ♦ 
0CNBA(12) * 
CYR (12) * 


X4* X5* XMRP* 

Y4* YS* YMRP* 

Z4* Z5* ZMRP* 
V(12) * 
CLBOa) * 
0CN9F(12) * 
CLR(12) * 

QQ(12) * 


DAMAX (12) *ORMAX(12) * 
CDA* CDR 


rLK 0010 

BLK 0020 
BLK 0030 
BLK 0040 
— BLK 0050 
BLK 0060 
IN THE BLK 0070 
COMPUTEBLK 0080 
BLK 0090 
BLK 0100 
BLK 0110 
BLK 0120 
BLK 0130 
« « « « plk 0140 
» » * » BLK 0150 
BLK 0160 
BLK 0170 
BLK 0180 
BLK 0190 
BLK 0200 
BLK 0210 
BLK 0220 
BLK 0230 
BLK 0240 
VY<12)* BLK 0250 
CNB(12)* BLK 0260 
CYA(12) * 8LK 0270 
CNR (12) BLK 0280 
BLK 0290 
BLK 0300 
BLK 0310 
BLK 0320 


DATA 

M*NS«KMAX 

•EPS0*MPT / 7* 3* 3* 0, 

00001* 12 / 


BLK 

BLK 

0330 

0340 

data 

TF / 

25.0 * 40.0 * 50.0 * 

60.0 * 65.0 * 70.0 

* 

BLK 

0350 

1 

DATA 

DATA 

DATA 

DATA 

X1,Y1*Z1 

X2*Y2*Z2 

X3,Y3*Z3 

X4»Y4*Z4 

75.0 ♦ 80.0 9 90.0 * 

/ 0.* 0. * -9*34 / 

/ 0.* 1.346* -6.68 / 

/ 0.* -1.346* -6.68 / 

/ 0.* 0. * 0. / 

100*0 * 110.0 « 140*0 

/ 

BLK 

8LK 

BLK 

BLK 

BLK 

BLK 

0360 

0370 

0380 

0390 

0400 

0410 

data 

X5*V5*Z5 

/O.* 0* • 0. / 



BLK 

BLK 

0420 

0430 

DATA 

XCG / 

23.345* 23.42 * 23.47 * 

23.52 * 23.545* 23.57 

* 

BLK 

0440 

1 

24.13 * 24.18 * 24.33 » 

24.535* 24.74 * 25.62 

/ 

BLK 

BLK 

0450 

0460 

DATA 

ZCG / 

-1.58 *-1.5847 *-1.59 14 *-1.5953*-!. 5979*-!. 60 

♦ 

BLK 

0470 

1 

-1 .4626*-! .455 *-1.440 *-l .4327*-! .4255*-! .400 

/ 

BLK 

BLK 

0480 

0490 

DATA 

XMRP*YMRP*2MRP / 21.6* 0.* -1.47 

/ 


BLK 

BLK 

0500 

0510 

DATA 

Q / 

.482E+4 * .987E+4 * 

.134E+5 * .174E^5 

♦ 

BLK 

0520 

1 


.194E+5 • *212E'»-5 * 

.226E4 5 * .23,3E>5 

9 

BLK 

0530 

2 


.217E+5 * .165E+5 * 

.117E+5 * .231E+4 

/ 

BLK 

BLK 

0540 

0550 

DATA 

S * 8REF 

/ 317.73 * 28.322 / 



BLK 

BLK 

0560 

0570 

DATA 

V / 

95.4 * 150. * 190. * 

241. * 272. * 305. 

* 

BLK 

0580 

1 


343. * 365. * 486. * 

612. * 768. t 1520* 

/ 

BLK 

0590 
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C 


C 


C 

c 


c 


c 


c 


c 


c 


c 


c 

c 


c 


c 


c 


c 


c 


c 


hata 

VY 

/ 

2. 

t 9. 

. 15. 

1 



AO, 

« 44. 

* 30. 

DATA 

F 

/ 

l.65F^6 * 1.76F+6 

1 



1.92F^<S * 1.94E^ft 

2 



2.02BE-^6 * 2. 

04E + 6 

OATA 

FS»M 

/ 

12*0. 

/ 


DATA 

cyb 

/ 

-1.66 

-1.68 

-1.70 

1 



-1 .97 

-1.92 

-1.93 

DATA 

CLR 

/ 

-.283 

-.285 

-.286 

1 



-.384 

-.356 

-.299 

DATA 

CNR 

/ 

.30? 

.315 

.325 

1 



.344 

.266 

.238 

DATA 

DCYR 

/ 

-.011 

-.012 

-.013 

1 



-.0165 

-.014 

-.0105 

DATA 

OCLR 

/ 

.0031 

.0032 

.0033 

1 



.0042 

.0035 

.0027 

DATA 

DCNBA 

/ 

.0064 

.0067 

.0074 

1 



.01 

.0088 

.0075 

DATA 

OCNBF 

/ 

-.004 

-.0044 

-.0046 

1 



-.0058 

-.005 

-.0044 

DATA 

CYa 

/ 

12*0. 

/ 


DATA 

CLA 

/ 

-.0430 

-.0458 

-.0487 

1 



-.0544 

-.0458 

-.0286 

DATA 

cna 

/ 

.0458 

. 0444 

.043 

1 



.0258 

.0244 

.0172 

DATA 

CYR 

/ 

.504 

• 408 

.462 

1 . 



.292 

.217 

.132 

DATA 

CLP 

/ 

.273 

.265 

.259 

1 



.206 

.186 

.105 

DATA 

CNR 

/ 

-.510 

-.489 

-.473 

1 



-.340 

-.254 

-.137 

DATA 

oamax 

/ 

40. 

40. 

40. 

1 



9.47 

8.91 

10.69 

DATA 

damax 

/ 

15. 

15. 

* 15. 

1 



9.47 

8.91 

10.69 

DATA 

DRMAX 

/ 

30. 

30. 

23.5 

1 



5.54 

5.23 

6.27 

OATA 

QQ 

/ 

•482E+4 • .987E+4 

1 



•194E+5 t .212E+5 

2 



.217F*5 , .165E*5 






8LK 

0600 

24. 

9 29. 

9 34 , 

9 

RLK 

0610 

0. 

9 0 , 

9 0. 

/ 

6LK 

0620 





RLK 

0630 

1.825F^6 f 1.1 

385E+6 

9 

8LK 

0640 

1.97E+6 t 1.98E^6 

• 

RLK 

0650 

2.06E+6 * 2. 

07E + 6 

/ 

RLK 

0660 





RLK 

0670 





RLK 

0680 





RLK 

0690 

-1.83 

9-1.99 

-2.05 

9 

BLK 

0700 

-2.03 

9-1.98 

-1.60 

/ 

BLK 

0710 





BLK 

0720 

-.291 

9-. 298 

-.326 

9 

BLK 

0730 

-.246 

9-. 196 

-.122 

/ 

RLK 

0740 





BLK 

0750 

.404 

9 .468 

.460 

9 

BLK 

0760 

• 269 

9 .207 

-.0284 

/ 

BLK 

0770 





RLK 

0780 

-.015 

9-. 016 

-.017 

9 

BLK 

0790 

-.008 

9-. 006 

-.004 

/ 

BLK 

0600 





BLK 

0810 

.0036 

9 . 0038 

.0042 

9 

BLK 

0820 

.0017 

9 .0014 

.001 

/ 

BLK 

0830 





BLK 

0840 

.0085 

9 .0094 

.0104 

9 

RLK 

0850 

.005 

9 .004 

• 0028 

/ 

BLK 

0860 





RLK 

0870 

-.0056 

9 - .006 

-.006 

9 

BLK 

0880 

-.0028 

9-. 0022 

-.0015 

/ 

BLK 

0890 





RLK 

0900 




• 

RLK 

0910 





RLK 

0920 

-r.0544 

9 -.0630 

-.0630 

9 

RLK 

0930 

-.0215 

9-. 0158 

-.00859 

/ 

RLK 

0940 





BLK 

0950 

.0358 

9 .0344 

.0301 

9 

BLK 

0960 

.00286 

9-.00286 

-.0114 

/ 

BLK 

0970 





BLK 

0980 

.394 

9 .319 

.300 

9 

BLK 

0990 

.0961 

9 .0749 

.0573 

/ 

BLK 

1000 





BLK 

1010 

.215 

9 .181 

.173 

9 

BLK 

1020 

.055 

9 .0406 

.0286 

/ 

BLK 

1030 





BLK 

1040 

-.388 

9-. 310 

-.345 

9 

BLK 

1050 

-.105 

9-. 077 , 

-.061 

/ 

RLK 

1060 





RLK 

1070 

40. 

9 25.1 1 

14.1 

9 

BLK 

1080' 

17.5 

9 33.64 1 

40. 

/ 

BLK 

1090 





BLK 

1100 

9 15* 

9 15. 

9 14.1 

.9 

BLK 

1110 

15. 

9 15. i 

15. 

/ 

BLK 

1120 





RLK 

1130 

.14.7 

9 8.19 1 

8.19 

9 

BLK 

1140 

10.23 

* 19.67 1 

30.0 

/ 

RLK 

1150 





BLK 

1160 

.134E+5 * .174E45 

9 

BLK 

1170 

.226E^5 • .233E^5 

9 

BLK 

1180 

.117E+5 9 .231E^4 

/ 

BLK 

1190 
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c 

DATA D«AXt SAt SR» COA* COR / 30. • 4*0. / 

C 

F.NO 


RLK 1200 
RLK 1210 
RLK 1220 
RLK 1230 



129 


1 



c INT 0010 

C INT 0020 

C SUBROUTINE INPUT(IGRAD*EPS*STEP^ISOf ICASE) INT 0030 

C INT 0040 

C IMT 0050 

C INT 0060 

C PURPOSE SUBROUTINE USED TO READ IN AND PRINT OUT THE INPUT DATA, INT 0070 

C INT 008 0 

C INT 0090 

C INPUTS rCASF = NO. OF CURRENT CASE. INT 0100 

C INT 0110 

C INT 0120 

C OUTPUTS IGRAO r ORDER OF GRADIENT METHOD TO BE USED. INT 0130 

C fPS s CONVERGENCE BOUND. INT 0140 

C STEP = STEP SIZE IF 1ST ORDER GRADIENT METHOD USED. INT 0150 

C IGO = INDEX CONTROLTNG SEQUENCE OF CASES. INT 0160 

C INT 0170 

C SUBROUTINES CALLED NONE INT 0180 

C IMT 0190 

C INT 0200 

0 jfgj 0210 

C 0220 

C INT 0230 

SUBROUTINE INPUT ( IGRAD*EPS*STEP» 160* ICASF) INT 0240 

C INT 0250 

C INT 0260 

1000 F0RMAT(7?AI) INT 0270 

1010 format { IHI *9X*4HCASE* I3*8X*72A1 / 10X*4H ) INT 0280 

1020 format (//5X,30HCOMPUTATION CONTROL PARAMETERS) INT 0290 

1030 FORMATnS*5X»2E10.3) INT 0300 

1040 p-ORMAT (5X*30H *10X*6lHUSE 1ST ORDERINT 0310 

1 GRADIENT METHOD ITERATION STEP SIZE = *E10.3 / 50X*3H ) INT 0320 

1050 format (5X*30H *10X*3lHUSF 2ND ORDERINT 0330 

1 gradient METHOD / 50X,3h ) INT 0340 

1060 format (58X.48HUPPER ROUND USED IN THE CONVERGENCE CRITERION = * INT 0350 

1 Flo. 3) INT 0360 

1070 format (5X. 17HTRAJECTORY POINTS ) INT 0370 

1080 F0PMaT(12I5) INT 0380 

1090 format (5X* 17H *13X*12I5 ) INT 0390 

non foRMAT(//5X*26HSYSTEM dynamics PARAMETERS) INT 0400 

1110 foRMaT<2F10.0) INT 0410 

1120 >^0RMAT (5X *26H t9X* INT 0420 

1 18H YAW BIAS TORQUE =*F11.1 / INT 0430 

2 40X*18HROLL BIAS TORQUE =*F11.1 ) INT 0440 

1130 format (//5X*32HPERF0RMANCE CRITERION PARAMETERS) INT 0450 

1140 FORMAT(7F10.0) INT 0460 

1150 format (5X*32H * INT 0470 

1 13X*5HW11 =*F7.2. 15X*5HW21 =*F7.2 / INT 0480 

2 50X*5HW12 =*F7.2*15X*5HW?2 =*F7.2 / INT 0490 

3 35X*9HWFIGHTIN6*6X*5HW13 = *F7 . 2 ♦ 1 5X *5HW23 =*F7.2 / INT 0500 

4 35X,9HFACT0RS *6X*5HW14 = *F7 .2* 15X » 5HW24 =*F7.2 / INT 0510 

5 50X*5HW15 =*F7.2* 15X*5HW25 =*F7.2 / INT 0520 

6 50X*5HW16 =*F7.?*15X*5HW26 =*F7.2 / INT 0530 

7 50X*5HW17 =*F7.2* 15X*5HW?7 =*F7.2 ) INT 0540 

1160 FORMAT (ID INT 0550 

1170 FORMAT <//5X,66H'» « * ERROR IN THE INPUT DATA — COMPUTER RUN TEINT 0560- 

IRMINATEO * * « ) INT q570 

C INT 0580 

dimension IDP(50) * ICP(50) * LINE(72) INT 0590 
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COMMON /CON/ M ♦ NS • KMAX « EPSO ♦ MPT 
COMMON /TRAJ/ JPT(12> * TF(1^) 

COMMON /SYST/ YBT? RRT 
COMMON /PEWF/ Wl(7>, W?(7) 


ENTER CASE IDENTIFICATION TITLE 
REAO(5*1000) (LINE(I) tl=l*72) 
write (6-1010) ICASE • (LINE( I) *1=1*72) 

*** ENTER COMPUTATIONAL CONTROL PARAMETERS 
WRITE (6*1020) 

RE80(S,1030) IGMAO * EPS * STEP 
lF(I(^RAO-2> 10*20*10 
10 WRITE (6* 1040) STEP 
OO TO 30 
20 WRITE (6*10S0) 

30 WRTTF(6*1060) EPS 

ENTER POINTS ALONG TRAJECTORY FOR COMPUTING TRIM 
WRITE (6*1070) 

READ(5*1080) ( JPT ( I ) * 1 = 1 *MPT) 

WRITE (6*1090) (JPTd) *I = 1*MPT) 

ENTER SYSTEM DYNAMICS PARAMETERS 
WRITE (6*1100) 

REAO(5*U10) YPT * RBT 
WRITE (6*1120) YRT * RRT 

*♦4 ENTER PERFORMANCE CRITERION PARAMETERS 
WR1TE(6*1130) 

REAO(5*1140> (Wl(I)*Ial*M) 

READ(5*1140) ' (W2(I) *I=1*M) 

WRITE (6*1150) Wl(l) * W2(l) * Wl(2) * W2(2) * Wl(3) ♦ W2<3) * 
1 Wl(4) * W2(4) t Wl(5) > W2(5) * Wl (6) * W2(6) » 

? Wl<7) 9 W2(7) 

ENTER END OF CASE CARO 
READ(S*1160) IGO 
TP(IGO-l) 200*220*200 
200 IF(IGO-2) 210*220*210 
210 WRITE (6*1170) 

CALL EXIT 
220 CONTINUE 
C 

RETURN 

END 


INT 0600 
INT 0610 
INT 0620 
INT 0630 
INT 0640 
INT 0650 
TNT 0660 
INT 0670 
INT 0680 
INT 0690 
INT 0700 
INT 0710 
INT 0720 
INT 0730 
INT 0740 
INT 0750 
INT 0760 
TNT 0770 
INT 0780 
INT 0790 
INT 0800 
INT 0810 
INT 0820 
INT 0830 
INT 0840 
INT 0850 
INT 0860 
INT 0870 
INT 0880 
INT 0890 
INT 0900 
INT 0910 
INT 0920 
TNT 0930 
INT 0940 
INT 0950 
INT 0960 
INT 0970 
INT 0980 
INT 0990 
INT 1000 
INT 1010 
INT 1020 
INT 1030 
INT 1040 
INT 1050 
INT 1060 
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INPUT Flow Dlogrann 



C * * — OUT 0010 

C OUT 0020 

C SUBHOUTINE OUTPUT < K *Lf M* T IME * DELTA »MPT > OUT 0030 

C out 0040 

C — — — —OUT 0050 

C . OUT 0060 

C PURPOSE SUBROUTINE USED TO PRINT OUT THE RESULTS OF THE PROGRAM. OUT 0070 

C , OUT 0080 

C OUT 0090 

G INPUTS K s NO. OF ITERATIONS. OUT 0100 

C L s NO. OF THE trajectory POINT. OUT OllO 

C M a NO. OF TRIM ANGLES. OUT 0120 

C TIME = FLIGHT TIME OF THE TRAJECTORY POINT. OUT 0130 

C delta = VECTOR OF TRIM ANGLES. OUT 0140 

C MPT = INDEX USED TO DETERMINE LAST TRAJECTORY POINT. OUT 0X50 

C OUT 0160 

C OUTPUTS NONE OUT 0170 

C — OUT 0180 

C OUT 0190 

C SUBROUTINES CALLED NONE OUT 0200 

C OUT 0210 

C - OUT 0220 

C *«*-»*************♦***♦♦*******♦#* ouj Q230 

Q «•»««««««««*<»«**««<»««**««»»«*«»««« o(jy 0240 

C OUT 0250 

SUBROUTINE OUTPUT (K*LtM*TlME*OELTAfMPT> OUT 0260 

C OUT 0270 

C OUT 0280 

1010 FORMAT (////5X*22HTRIM. deflection ANGLES * OUT 0290 

1 /5X*22H *7X»13HTRAJ. FLIGHTt26Xt OUT 0300 

2 5HDELTA*32X*6HN0. OF / 35X*85HPT. TIME (1) (2) (3) OUT 0310 

3 (4) (5) (6) (7) ITERATIONS ) OUT 0320 

C1020 FORMAT(35X*I3*2X*F6.1tlX*7F0.2*7X*I5) OUT 0330 

1020 FORMAT(3lX*lH.*3X*l3*2XfF6.1f lX*7F8.2t4X*lH.t2X*I5) OUT 0340 

1030 FORMAT! 48X*55H TOP YAm/ PITCH YAW PITCH AILERON RUOUT 0350 

lODER * /48X»40H<— — OR0ITER — — — SRM > ) OUT 0360 

1040 F0RMATOlX«77H .••••«... OUT 0370 

) OUT 0380 

1850 FORMAT OIXrlH.fTSX.lH. > OUT 0390 

C OUT 0400 

DIMENSION DELTAU) * ANGLE(lO) OUT 0410 

DATA HAD / 57.2957795 / OUT 0415 

C OUT 0420 

C OUT 0430 

TF(L-I) 20, 10*20 OUT 0440 

10 WRITE(6,1010) out 0450 

WRrTE(6,1040> OUT 0460 

WRITE(6f 1050) OUT 0470 



20 IF(K) 40*30»30 

30 OP 36 1=1 »M 

36 ANGLE<1) = RAO » OELTAU) 

WRITE(6* 1020) L * TIMt » ( ANGLE ( I ) * 1 = 1 *M) » K 

40 IF(L-MPT) 60^60»60 
50 WRITE(6»1060) 

WRITE(6»1040) 

WRITE(6*1030) 

60 RETURN 
END 


OUT 0480 
OUT 0485 
OUT 0486 
OUT 0490 
OUT 0500 
OUT 0510 
OUT 0520 
OUT 0530 
OUT 0540 
OUT 0550 





oo ooo on on 


SUrtrtOUTlNt GWa01(»v*LiT I Mt * OELT A ♦ L AMDA ♦ I GRAD* EPS * STEP ) 


PURPOSE SUdPOUTlNt FOR COMPUTING THE DEFLECTION ANGLES USING THE ONE 


inputs 


OUTPUTS 


1ST OPDEH GRADIENT METHOD, 

K = NO. OF iterations. 

L = NO. OF the trajectory POINT. 

TIME = FLIGHT TIME OF THE TRAJECTORY POINT. 
delta = INITIAL guess" OF TRIM ANGLES. 

LAMOA = INITIAL GUESS OF LAGRANGE MULTPLIERS. 
IGkAU = i 

EPS = CONVERGENCE HOUND. 

STEP = STtH SIZE. 

delta = VECTOR OF TRIM ANGLES. 

LAMDA = VECTOR OF LAGRANGE MULTIPLIERS. 


SUHROUTINES CALLED SYSTEM * COST ♦ MCPY ♦ CCUT * MINV * 6MPRD .ONE 

ONE 

ONE 

««.*«-i^««tto-*-»«»*«««-*^*«’**'*********** ONE 
»»»#»*««»»* *»<*^*”*^***************** ONE 

ONE 

subroutine GRADl (»\tL*TIME*l)ELTA*LAMDA*IGRAO*EPS*ST£P) ’ ONE 

ONE 

ONE 

1000 F0RMAT.(//SX.TSH<»* warning 1ST ORDER GRADIENT ALGORITHM USED THONE 

IE MAX. NO, OF ITERATIONS* 14 /i20X*6HNORM = * E 1 0 .3 * 1 OXtSHEPS =*E10.3)ONE 

ONE 

. * TYPE AND STORAGE ALLOCATION •••.•••• .ONE 

HEAL LAMOA * NORM 

DIMENSION DELTA(IO)* LAMDA(b)» BX(60)* BU(60>* RX(10)» RUdOl* ONE 

1. XnO) * DX(IO) * OU(IO) * LBUO) * MB(IO) ONE 

COMMON /ARRAY/ AV(6)* BV(6>* 6M(60)* BT(6*60)* RS* RV(XO)* RM(100)ONE 
COMMON /CON/ M * NS * KMAX ♦ EPSO * MPT ONE 

COMMON /TRAJ/ JPlilk) * TF(12) 

ONE 

ONE 

TEST, whether THIS trajectory POINT IS TO BE USED ONE 

IF(JPT(D) 5*l#b 

1 K = -1 ONE 

GO TO 130 , 

ONE 

COMPUTE the time OF FLIGHT , ONE 
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oo oo oo no on onoo n no 


Te (L) 


b TlMt = 


st&ht Initial itlhation 

K = 0 

MNS = M « NS 


10 CONTINUt 

*** COMPUTE GPAOIENT TERMS CORRESPONDING TO SYSTEM DYNAMICS 
CALL SYSTEM (KfLfNSfMtOELTA* IGRAD) 

PAhTITION The MATRIX BM INTO MATRICES BX AND BU 
IF(MNS) 30*30*40 
30 CALL MCPY {BM*BX*NS*M*0) 

60 TO bO 
40 J = NS + 1 

CALL CCUT (8M* J*BX*BU*NS*M*0) 

*** COMPUTE THE INVERSE OF THE MAThIX BX 
SO CALL MINV (BX*NS*OET*LB*Me) 


*** COMPUTE VECTOR X 
00 SO I=1*NS 

60 DU(1) = - AV(I) - BV(I) 

CALL GMPRU (BX*DU9DX«NS*NS* 1 ) 
DO 6b I=I*NS 
6b X(I) = X(I) ♦ UX(I) 


*** COMPUTE GRADIENT TERMS CORRESPONDING TO PERFORMANCE CRITERION 
00 70 I=I*NS 
70 delta (1) = X(I) 

IF(MNS) 130*130*80 
80 CALL COST{K*L*M*OELTA*I6RAO) 

44* PARTITION THE vector RV INTO VECTORS RX AND RU 
J = NS ♦ I 

CALL CCUT (RV* J*RX*RU*1*M«0) 

*4* COMPUTE the VECTOR LAMOA 

CALL GMPRD(RX*BX*LAMDA*1#NS*NS> 

*4* COMPUTE The new estimate of delta 

CALL. GMPRO (LAMDA*BU*0U*1*NS*MNS) 

NORM =0. 

DO 90 I=1*MNS 

OU(I) = ( DU(I) - RUil) ) 4 STEP 
NORM = NORM ♦DU <1)442 
90 delta (NS+I) = DELTA(NS^I) ♦ OU ( I ) 


ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 


0490 

OSOO 

OblO 

0S20 

0S30 

0540 

0550 

0560 

0570 

0580 

0590 

0600 

0610 

0620 

0630 

0640 

0650 

0660 

0670 

0680 

0690 

0700 

0710 

0720 

0730 

0731 

0732 
0740 
0750 
0760 
0770 
0780 
0790 
0600 
0810 
0820 
0830 
0840 
0850 
0860 
0870 
QQSO 
0890 
0900 
0910 
0920 
0930 
0940 
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DO on no 


00 9b I=1*NS 

9b N0^^M = NuKm + UA{I)**2 

TEST IF THE NEW ESTIHATtS ARE SUFFICIENTLY ACCURATE 
IF<N0RH°£RS) 13O»l309l00 

«•«* CHECK FOR EXCESSIVE NUiWbfER OF ITERATIONS 
100 IF{K-KMAA) 1 1 Of 1<^0» lifO 

perform another iteration 
110 K = K + 1 
GO TO 10 

li?0 WRITE (fet 1000) Kf NORM? EPS 
IJO RETURN 
ENO 


ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 

ONE 


0941 

0942 
U9b0 
0960 
0970 
0960 
0990 
1000 
1010 
1020 
1030 
1040 
1050 
1060 
1070 
1080 
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GRADl Flow Diagram 







GRAD 1 Flow Diagram (Confinued) 
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C QQJQ 

^ TWO 0020 

C SU 8 ^>OUTIiMe (K,L, I XMtf OELTAf LAMOAt IGHADfEPS) TWO 0030 

^ TWO 0040 

C ^ TWO 0050 

C TWO 0060 

C PU^^POSE SUH-iOUTINE COMPUTING THE DEFLECTION ANGLES USING THE TWO 0070 

C ^^JO OP[)ER GRADIENT METHOD. TWO 0080 

C TWO 0090 

C INPUTS K = NO. OF iterations. . TWO 0100 

C L = NO. OF THE TRAJECTORY POINT. TWO OUO 

c time = flight TIME OF THE TRAJECTORY POINT. TWO 0120 

C DELTA = INITIAL GUESS OF TRIM ANGLES. TWO 0130 

C LAMDA = INITIAL GUESS OF LAGRANGE MULTPLIERS. TWO 0140 

C IGRAf.) = 2 TWO 0150 

C EPS = convergence bound. two 0160 

^ TWO 0170 

C OUTPUTS DELTA = VECTOR OF TRIM ANGLES. TWO 0180 

C LAMDA = VECTOR OF LAGRANGE MULTIPLIERS. TWO 0X90 

C two 0200 

C SUBROUTINES CALLED SYSTEM * COST , SINV « MXOUT * MPRO f TWO 0210 

C GmTRA « GMPRD t GMSYMM 9 MSTR • GMSU0 .. TWO 0220 

^ TWO 0230 

C * TWO 0240 

C *«»*•»*'»*<>'»*«««»**•»*«*«**♦♦*♦#**** TWO 0250 

^ TWO 0260 

SUBROUTINE GRA02(K*L*TIME90ELTA»LAHOA*I6RADtEPS) TWO 0270 

TWO 0280 
TWO 0290 

1000 FORMAT (//bXf5sH*« ERROR «« MATRIX R IS NOT POSITIVE DEFINITE TWO 0300 

1 ♦ 13 * bX *5HEPS =9E12.3 /) TWO 0310 ' 

1010 format (//5X965H#* WARNING ** LOSS OF SIGNIFICANCE IN INVERTING HATWO 0320 
ITRIX R K =9I3«5X95HEPS =*E12.3 /) TWO 0330 

1020 format (//5Xt57Hfr<» ERROR MATRIX 0R8 IS NOT POSITIVE DEFINITE TWO 0340 

1 K = 9 1 3« bX 9 bHEPS ~9£12.3 /) TWO 0350 

1030 format (//5X967h«* WARNING ** LOSS OF SIGNIFICANCE IN INVERtiNG MATWQ 0360 

. ITRIX BR8 K =f I3*bX,5HEPS =.E12.3 /) TWO 0370 

1040 FORMAT (//bX968H«* WARNING ** 2ND ORDER GRADIENT METHOD USED MAX. TWO 0380 

INO. OF lTERATI0NS9l395Xt5HEPS =9E12.5*5X96HN0RM =9E12*5 /) TWO 0390 

lObO FORMAT (/10X99HMATRIX R ) TWO 0400 

1060 FORMAT (/10X.20HMATRIX R (INVERSE)) T«0 0410 

1070 FORMAT (/10X99HMATRIX 8 ) TWO 0420 

1080 format (/I0X9 1 InmaTRIX 8R8 ) XWO 0430 

C TWO 0440 

REAL LAHOA 9 LAM 9 NORM TWO 04S() - 

DIMENSION OELTA(10)9 LAMUA(6)* OELdO)* LAM(6)» HO (1 0 ) * HLU 0) » TWO 0460 

1 R(IOO)* 0(60)9 BR(60)* BRB ( 36 ) 9 D (60 ) 9 X(10)» Y(10) TWO 0470 

COMMON /ARRAY/ AV(6)t BV(6>9 BM(60>9 BT(6*60)t RSt RV ( 1 0) t RM ( 1 00 ) TWO 0480 



COMMON /CON/ M » NS ♦ KMftX ♦ tPSO ♦ MPT 
common /THaJ/ JPT(12) ♦ TFMie) 

FOOIVALFNCF. (B(l)*t<M(l)) 

c 

c 

Q test whETHEk this thajectopy point is to be used 

IF { JPT (L) ) b» 1 f S 
1 K = ,~J 
GO TO IGO 
C 

Q »■»* cOMPUTt THc time of flight 
S time = TF(L) 
c 

K' = 0 


. c 

c- - 

C OPTION FOH OTSPEGAhOING AILEPON 

IF(JPT(L)+ii) 
b IGPAO = - IGPAD 

'■ C 

1 1; CALL S Y S T E M ( K < L * NS « M ♦ DEL T A f 1 GP AO ) 
call cost (K »L*m*OELTA» IGPAD) 

c ' 

c COMPUTE The OEPIVATIVE of the HAMILTONIAN WITH RESPECT TO DELTA 

DO I = lf^ 

HO ( I ) = P^/ ( I ) 

00 ^0 J=1^ ^IS 
ji = j + n-i)">NS 

i?0 HU { I) = HO ( I ) + LAMUA ( J) ( JI ) 

c • 

c compute the DERIVATIVE OF THE HAMILTONIAN WITH RESPECT TO LAMDA 

00 JO J=1*nS 

30 HL(J) = AV(J) + BV(J) 

C 

c COMPUTE Thl ;^N0 derivative of the HAMILTONIAN R = HOD 
M3 = M^^(M+l)/^ 

DO ^0 1=1 
R t I ) = PM U ) 

DO 40 J=i f NS 

40 R(I> = PII) + L AmOA ( J) *HT ( J* I) 

C 

C COMPUTE THE ?N0 DERIVATIVE OF THE HAMILTONIAN B = HLO 

c • 

c (SEE EOUIVALENCE STATEMENT ) 

C 

c COMPUTE INVERSE OF MAThIX R 

CALL SINV (R*M»EPSU* IEk) 


TWO 

0490 

TWO 

0500 

TWO 

0510 

TWO 

0520 

TWO 

0530 

TWO 

0540 

TWO 

0550 

TWO 

0560 

TWO 

0570 

TWO 

0580 

TWO 

0590 

TWO 

0600 

TWO 

0610 

TWO 

0620 

TWO*0621 

- TWO*0622 

TWO*0623 

TWO*0624 

TWO*0625 

- TW0*Q626 

TWO*0627 

TWO 

0630 

TWO 

0640 

TWO 

0650 

TWO 

0660 

TWO 

0670 

TWO 

0600 

TWO 

0690 

TWO 

0700 

TWO 

0710 

TWO 

0720 

TWO 

0730 

TWO 

0740 

TWO 

0750 

TWO 

0760 

TWO 

0770 

TWO 

0780 

TWO 

0790 

TWO 

080Q 

TWO 

0810 

TWO 

0820 

TWO 

0830 

TWO 

0840 

TWO 

0850 

TWO 

0860 

TWO 

0870 

TWO 

0880 

TWO 

0890 
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IF(IEr^) 60*7 0*b0 
50 Wf^iTF (b* 1000) K 9 EPSU 
WRITE (b* 1 ObO ) 

CALL .^AUU r ( I *K9M9f«l9 1 too * 132< 1 ) 
CALL EXIT 

60 WRrTE(6tlOiO) K t tPSO 
WRITE(6tl050) 

CALL HXOUT (l9R*MtMtl*60*132tl) 

compute matrix hr 

70 CALL MPNO (BtH*HrttNSfM*Ot 1 tM) 


*«» compute matrix BRb 
CALL GMTRA (HtOtNS*M) 

CALL (iMPRU ( HRtOtHRrt f NbtMfNS ) 

call GMSrMM (BRBtDtNS) 
call msth ( d t brb » ns ♦ 0 t 1 ) 

*** compute inverse of matrix brb 
call SINV (BRBtNStEPSOt IER) 

IF(IER) 809100*90 
60 WRI TE (6* 1020) K 9 EPSU 
WRITE (691060) 

CALL HXOUr(l9R9M9M9l9609l329i) 
WRITE (6*1070) 

CALL MXOUT ( 1 9fl9NS9M90960* 132* 1 ) 
WRITE (6*1060) 

CALL MXl)(JT ( 1 9 HRB 9 NS 9 NS 9 1 *60 *132*1) 
CALL EXIT 

90 WRITE (6* 1030 > K * EPSO 
WRITE (6*1060) 

CALL MXOUT ( I 9 R 9 M 9 M 9 1 *60* 132* 1 ) 
WRITE (6*1070) 

CALL MXOUr (l9B9NS*M*U9609l32*l) 
WRITE (6*1080) 

CALL. MXOUT ( 1 *HR8*NS*NS* 1 *60 * 132* 1 ) 


option FOR OISREGARUIN6 1ST TkIM EOUALITY CONSTRAINT 
EQUATION RE(;)UIHING ZERO NET FORCE IN Y-OIRECTION 
lOQ IF(JPT(L)> 95*96*96 

95 AV(1) = HRB (2) «AV (2). ♦ BRB(4)«AV(3) 

HL(U = HV(l) - AV(1)/HK8(1) 

96 continue 


■»*» COMPUTE CORRECTION TO LAMOA 


TWO 

0900 

TWO 

0910 

TWO 

0920 

TWO 

0930 

TWO 

0940 

TWO 

0950 

TWO 

0960 

TWO 

0970 

TWO 

0960 

TWO 

0990 

TWO 

1000 

TWO 

1010 

TWO 

1020. 

TWO 

1030 

TWO 

1040 

TWO 

1050 

TWO 

1060 

TWO 

1070 

TWO 

1060 

TWO 

1090 

TWO 

1100 

TWO 

1110 

TWO 

1120 

TWO 

1130 

TWO 

1140 

TWO 

1150. 

TWO 

1160 

TWO 

1170 

TWO 

1180 

TWO 

1190 

TWO 

1200 

TWO 

1210 

TWO 

1220 

TWO 

1230 

TWO 

1240 

TWO 

1250 

TWO 

1260 

TWO 

1270 

- TW0*1271 
TW0^1272 
TW0*1273 
TW0*I274 
TWO* 1275 
TWO* 1276 
TW0*1278 

- TWO* 1278 
TW0*1279 

TWO 

1280 



nr) ' on DO 


call bMPHD (BRfHD»X»NS»M* 1 ) 

CALL wPkO ( HBB ♦ X f Y f NS f NS • 1 • 0 ♦ 1 ) 
CALL i-lPHlO ( BBH ♦ HL f X ^ NS* NS » 1 T 0 < 1 ) 
CALL (Xf Y»LAM*NS* 1) 

»««• COMPUTE COPPfcCTION TO DELTA 
call GMPKO ( Y^BBf DEL* 1 ♦Nb^M) 

CALL MPHO (P *HD*Y»M»Mt 1 < 0* 1 ) 

CALL (jMSUB (DLL * Y ♦DEL*M* 1 ) 

CALL 6i«PHU(X»BR»Y^l»NSfM) 

CALL GMSUb (l>EL* Y»DEL*M* i ) 

COMPUTE NE# estimate OF DELTA 
NOPM = 0* 

00 TlO I=1»M 

N0B(4 = NORM ♦ 0EL(1)^‘^*^ 

1 1 0 DELt A (I) = DEL T A ( I) + i?EL ( I ) 

compute NEW ESTIMATE OE' LAMDA 
DO 120 J=1*NS 
NOb4 = NOPM + LAM(J)*»2 
12U LAMOA(J) = LAMUA(J) ♦ LAM(J) 
IF{^ORm-EPS) 160»lbO*lJ0 
130 IF(K-t<.MAX) 140*lSO*lb0 

140 K = :k + 1 
60 IjO 10 

150 WRITE (6* 1040 ) K * FPSO » NORM 
160 RFTURN 
END ; 


TWO 1290 
TWO 1300 
TWO 1310 
TWO 1320 
TWO 1330 
TWO 1340 
TWO 1350 
TWO 1360 
TWO 1370 
TWO 1380 
TWO 1390 
two’ 1400 Y 
TWO 1410 
TWO 1420 
TWO 1430 
TWO 1440 
TWO 1450 
TWO 1460 
TWO 1470 
TWO 1480 
TWO 1490 
TWO 1500 
TWO 1510 
TWO 1520 
TWO 1530 
TWO 1540 
TWO 1550 
TWO 1560 
TWO 157Q 
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GRAD 2 Flow Diagram (Continued) 
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on on 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

100 


SYS 0010 

SYS 0020 

SUBHOUTIi'^E SYSTEM (K*L*NS*M*D£LTA* IbrtAO) SYS 0030 

SYS 0040 

. SYS 0050 

SYS 0060 

PUh{POSK SUHPOUriNt FOK COMPUTING ThE COEFFICIENTS IN THE SYS 0070 

EUOATIONS OF THE LATEHAL OYNAMICS DEFINING THIN* SYS 0080 

ALSO evaluates THE CORRESPONDING DERIVATIVES REQUIRED SYS 0090 
«Y The gradient METHODS. SYS OlOO 

SYS 0110 

inputs K = NO. OF iterations. SYS 0120 

L = NO. OF THE trajectory POINT. SYS 0130 

NS = NO. OF TRIM EQUATIONS. SYS 0140 

M = NO. UF trim angles. SYS 0150 

delta = VECTOR OF TRIM ANGLES. SYS 0160 

IGHAD = ORDER OF GRADIENT METHOD TO 8E USED. SYS 0170 

SYS 0180 

SUBROUTINES CALLED NONE SYS 0190 

SYS 0200 

SYS 0210 

•tt'-tt**»'a-*'M^«*<j’tt*<><J'»*<>******^****'^*** SYS 0220 
•M’*'tt**Oil-‘»'tt«-*'il-'i»'**‘tti*’***^*********'*^** SYS 0230 

SYS 0240 

SUBROUTINE SYSTEM (K*L*NS*M*0ELTA* IGRAO) SYS 0250 

SYS 0260 
SYS 0270 

DIMENSION delta U) ' SYS 0280 

COMMON /array/ AV(6)t bV(6)* BM(60)» BT(6*60)* RSf RV(10)t RM(100)SYS 0290 
COMMON /SYST/ YbT * R6T . XI* A2* X3* X4t X5* XMRP* SYS 0300 

1 S* HREF* Yit Y2. Y3* Y4* Y5* YMRP* SYS 0310 

^ Zlt 22* 23* 24* 25* 2MRP* SYS 0320 

3 XCG(12)* ZCG(12)* U(12)* V<12)* V Y ( 1 2 ) * SYS 0330 

4 F(12)* FSHM(12)* CY0(12)* CL8(12)* CN8<12)* SYS 0340 

5 0CYBU2). 0CLB(12)* DCNBA(12)* DCNBF(12)* CYA(12>* SYS 0350 

6 . CLA(12)t CNA(12). CYR(12)* CLR(12). CNR(12) SYS 0360 

DATA RAD / 57.2957795 / SYS 0365 

SYS 0370 
SYS 0380 

IF(K) 300*100*300 SYS 0390 

COMPUTE vector A SYS 0400 

CONTINUE SYS 0410 

QS = Q(L) S SYS 0420 

QSB = QS « BREF SYS 0430 

BETA = ARSIN(VY <L)/V (L) ) SYS 0440 

SYS 0450 

CYBCG = CYB(L) ♦ (UCYB (L) ) *«AD SYS 0460 

CLBCG = CL8(L) ♦ (OCLB (L) ) *RAD SYS 0470 
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CNBCG = CnB(L) ♦(l>CNtJA(L) ♦ OCNeF(L))*HAO 
C 

CLHClV = GLBC6 + C YBC<i* ( ZCG (L ) -ZMKP) /BREf 
• CNtiCG = CNBC6 - CYHGG« UC6 (L) -XMRP ) /B«EF 
C 

AV(1) = OS ^ CYbCG « beta 

AV{^) = USB « CLBCG * BETA ♦ RBT 

AV(3) = OSB « CNBCB « BETA ♦ YBT 

**« compute coefficients in vector B 

WAD = S7.ii95779S 
Cl = C0S(18,/WAD> 

51 = SIN(lo./RAD) 

C2 = C0S<12./WAD) 

52 = SIN(l2,/RAlj) 

C3 = C0S(3. S/RAD) 

53 = S1N(3.S/WA0) 

C4 = COS(lb./HAU) 

54 = SIN(ib./WAO) 

CLACG = CLA(L) ♦ C YA (L ) * ( ZC6 (L ) -2MRP) /8REF 
CNAC6 = CNA(L) - CYA(L^<^(XCG(L)-XMRP)/BR6F 

CLRC6 = CLW(U ♦ CYW(L)*(ZCG(L>-2MRP)/8REF 

CNRCG = CnW<U - CYR(L)<^(XC6(L)-XMRP) /BREF 

BM(i) = F (L) * CI 

BM(2) =-F(L) * Cl » (21 - ZCG(L>) 

RM(3) s F(L) * Cl * (XI - XC6(D) 

BMI4) = 2. » F(L) « C2 * C3 

BM(b> =-2. * F(L) « C2 C3 * (Z2 - ZCG(L>) 

8M(6) = 2. F(L) « ( (X2-XCG(L) )*C3 - Y2*S3) * C2 

BM(7) = 2. » F (L) ■» S2 * S3 

BM(8) = 2, * F(L> * (Y2*C2 - ( Z2-ZC6 (L ) ) *S2*S3) 

8M(9) = 2, * F(U « ((Y2*C3 ♦ (X2-XC6 (D ) *S3) * S2) 

BM(IQ) X 2. * FSRM(L) » C4 
BMUl) 5-2. ft FSRM(U * C4 * (Z4 - ZC6<D) 

8M(12) 5 2. * FSRM(L) ♦ ( ( X4-XCG ( L H *C4 - Y4*S4) 
BM(13) 5 0. 

8M(14) = 2. ft FSRM(L) » Y^ 

8M(ib) =0. 

8M(16) 5 as ♦ CYA(L) 

8M(17.) = QSB * CLACG 

BM (18) s QSB * CNACG 

8M(19) 5 OS ♦ CYH(L) 

8M(20) 5 USB * CLRCG 

BM(21) = QSB ♦ CNRCG 


SYS 0480 
SYS 0490 
SYS 0500 
SYS OSIO 
SYS 0S20 
SYS 0530 
SYS 0540 
SYS 0550 
SYS 0560 
SYS 0570 
SYS 0580 ‘ 
SYS 0590 
SYS 0600 
SYS 0610 
SYS 0620 
SYS 0630 
SYS 0640 
SYS 0650 
SYS 0660 
SYS 0670 
SYS 0680 
SYS 0690 
SYS 0700 
SYS 0710 
SYS 0720 
SYS 0730 
SYS 0740 
SYS 0750 
SYS 0760 
SYS 0770 
SYS 0780 
SYS 0790 
SYS 0800 
SYS 0810 
SYS 0620 
SYS 0830 
SYS 0840 
SYS 0850 
SYS 0860 
SYS 0870 
SYS 0880 
SYS 0890 ! 
SYS 0900 ' 
SYS 0910 ■ 
SYS 0920 
SYS 0930 1 
SYS 0940 
SYS 0950 
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««« OPTION FOP DISPtGAKUlNO AILERON 
IF(IGPAO) i?90*300«J0U 
2S0 = 0. 

HM(17) = 0, 

HM(IH) = 0* 

I6PAD = - IGPAD 


«*“» COMPUTE VECTOR b 
300 CONTINUE 

00 31U I=i»NS 
BV(I) = 0, 

DO 310 J=1»M 
IJ = I ♦ <J“1)*NS 

310 BV(I) = BV(I> ♦ BM ( 1 J) *DELTA ( J) 

COMPUTE The 1ST derivative of vector 8 

{ CONSTANT matrix COMPUTED ABOVE ) 


IF(IGRAD-2) B00*b00*b00 

COMPUTE THE 2ND DERIVATIVE OF VECTOR B 
BOO CONTINUE 

m 2 = M*(M+i)/2 
DO 510 J=ifNS 
00 510 1=1 *M2 
BIO BT ( I ) = 0. 

600 return 

END 


- SYS*0951 
SYS*0952 
SYS«0953 
SYS*0954 
SYS*0955 
SYS*0956 
SYS*t)957 
« SYS*0958 
SYS*0959 
SYS 0960 
SYS 0970 
SYS 0980 
SYS 0990 
SYS 1000 
SYS 1010 
SYS 1020 
SYS 1030 
SYS 1040 
SYS 1050 
SYS 1060 
SYS 1070 
SYS 1000 
SYS 1090 
SYS 1100 
, SYS 1110 
SYS 1120 
^ SYS 1130 
SYS lUO 
SYS 1150 
SYS 11$0 
SYS 1170 
SYS 1180 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CST 

CST 

SU^^ROUTINE COST (K*L»Mf DELTA* IGRAD) CST 

CST 

; CST 

CST 

PURPOSE SUBROUTINE FOR COMPUTING THE COEFFICIENTS IN THE CST 

performance CRITERION. CST 

ALSO EVALUATES THE CORRESPONDING DERIVATIVES REQUIRED CST 
BY THE GRADIENT METHODS. CST 

CST 

INPUTS K = NO. OF ITERATIONS. CST 

— L = NO. OF THE TRAJECTORY POINT. CST 

M = NO. OF TRIM angles. CST 

delta = VECTOR OE TRIM ANGLES. CST 

IGRAO = ORDER OF GRADIENT METHOD TO BE USED. CST 

CST 

SUBROUTINES CALLED NONE CST 

CST 

CST 


■tt»«»<^«»»<^««*<>«*********<^* **««»*** CST 

CST 

SUBROUTINE COST {K*L*M*0ELTA*I6RAD) CST 

CST 

CST 

DIMENSION C(7) * DELTAU) CST- 

COMMON /array/ AV(6)* BV(6)* BM(60)* BT(6*60)* RS* RV(10)* RM(100)CST 
COMMON /PERF/ Wl(7)* W2(7)* DAMaX (12) *ORMAX (12 ) ♦ QU2)* ’ CST 


1 DMAX* SA* SR* COA* COR CST 

CST 

CST 

IF <K) 100*200*300 CST 

CST 

COMPUTE MINIMUM VALUE OF THE PERFORMANCE CRITERION CST 

100 CONTINUE CST 

RS =0. CST 

DO no 1 = 1 *M GST 

110 RS = RS ♦ C(I) * 0ELTA(1)**2 CST 

RS = RS / 2. CST 

60 TO SOO CST 

CST 

**«• COMPUTE THE COEFFICIENTS IN THE PERFORMANCE CRITERION CST 

200 CONTINUE CST 

DO 210 1 = 1 *S CST. 

210 C(I) = (Wl (I)/DMAX)**2 ♦ W2 ( I ) **2 CST 

C(6) = (W1 (6) /DAMAX (L) ) «^«2 ♦ ( W2 ( 6 ) ( L ) *SA*CDA ) **2 CST 

C(7) = (Ml (7) /ORMAX (L> ) **2 ♦ ( W2 ( 7) *Q (D *SR*CDR) **2 CST 


i 

I 


OOlO 
0020 
0030 
0040 
0050 
0060 
0070 
0080 
0090 
0100 
0110 
0120 ’ 
0130 
0140 
0150 
0160 
0170 
0180 
0190 
0200 
0210 
0220 
0230 
0240 
0250 
0260 
0270 
0260 
0290 
0300 
03 io 
0320 
0330 
0340 
0350 
0360 
0370 
0380 
0390 
0400 
0410 
0420 
0430 
0440 
0450 
0460 
0470 
0480 


152 



n n n no 


CST 0461 
CST 0482 
CST 0483 
CST 0484 
CST 0490 
CST 0500 
CST 0510 
CST 0520 
CST 0530 
CST 0540 
CST 0550 


C *** COMPUTt THE 1ST DERIVATIVE OF Th£ PERFORMANCE CRITERION CST 0560 

300 continue CST 0570 

DO 310 I-lfM CST 0560 

310 RV{I) = C(I) * OELTA(I) CST 0590 

(j CST 0600 

IF(I6RA0-2) 500f400*500 CST 0610 

CST 0620 

**« COMPUTE THE 2ND DERIVATIVE OF THE PERFORMANCE CRITERION CST 0630 

400 continue CST 0640 

CST 0650 

- constant matrix computed ABOVE ) CST 0660 

CST 0670 

500 return GST 0660 

END CST 0690 


DO 212 I=1*M 
IF(wl(D) 2U*212»212 

211 CU) = !• / W1{I)**2 

212 continue 

M2 = /^ 

DO 220 1=1 *M2 
220 RM(I) =: 0. 

DO 230 1=1 »M 
II = I*( I^i)/2 
230 RMdl) = C(U 



C0ST Flow Diagram 
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/ 


SUBROUTINE GMSYMM ( A »8f N) 


PURPOSE CO^iPUTES A SYMMETRIC MATRIX B FROM A SQUARE MATRIX 
— ACCORDING TO 

B = ( A ♦ A*) / 2 

INPUTS A a SQUARE MATRIX (STORAGE MODE =0). 

QP COLS, IN A AND B, 

OUTPUTS B = SYMMETRIC MATRIX FORMED FROM A (STORAGE 

mode s 0), 

SUBROUTINES CALLED NONE 


»««««*«*»»«*««««««« 
subroutine GMSYMM(4»8*N> 


DIMENSION A(l) • R(l) 


N1 = N - 1 
IF(Nl) 20t20f5 
S DO 10 JalfNl 
J1 s J ♦ I 
DO 10 IsJltN 
IJ = (J-l)*N ♦ I 
JI = (I-l)*N ♦ J 
0(IJ) = 0,5 * (A(IJ) ♦ A{J1>) 
10 R(JI) = BdJ) 

20 DO 30 I«1*N 

TJ = ri-l)*N ♦ 1 
30 R(IJ) 3= A(IJ) 

RETURN 

END 



0010 

SYM 

0020 

SVM 

0030 

SYM 

0040 

•— — SYM 

0050 

SYM 

0060 

A SYM 

0070 

SYM 

0080 

SYM 

0090 

SYM 

0100 

SYM 

0110 

SVM 

0120 

SYM 

0130 

SYM 

0140 

SYM 

0150 

SYM 

0160 

SYM 

0170 

SYM 

0180 

SYM 

0190 

* * SYM 

0200 

* * SYM 

0210 

SYM 

0220 

SYM 

0230 

SYM 

0240 

SYM 

0250 

SYM 

0260 

SYM 

0270 

SYM 

0280 

SYM 

0290 

SYM 

0300 

SYM 

0310 

SVM 

0320 

SYM 

0330 

SYM 

0340 

SYM 

0350 

SYM 

0360 

SYM 

0370 

SYM 

0380 

SYM 

0390 

SYM 

0400 

SVM 

0410 

SYM 

0420 
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c 



MCPY 

10 

c 



.•••MCPY 

20 

c 



MCPY 

30 

c 

SUBRCUTINE J^CPY 


MCPY 

40 

c 



MCPY 

50 

c 

PURPOSE 


MCPY 

AO 

c 

COPY ENTIRE MATRIX 

. 

MCPY 

70 

c 



MCPY 

80 

c 

USAGE 


MCPY 

90 

c 

CALI MCPY (AtR»N»M*MS» 


MCPY 

100 

c 



MCPY 

110 

c 

CESCRIPTICN CF PARAMETERS 


MCPY 

120 , 

c 

A - NAME CF INPUT MATRIX 


MCPY 

130 

c 

R - NAME CF OUTPUT MATRIX 


MCPY 

140 

c 

N - NUMEER OF ROWS IN A CR R 


MCPY 

150 

c 

M - NUMBER OF COLUMNS IN A OR R 


MCPY 

160 

c 

MS - ONE DIGIT NUMBER FOR STORAGE 

MODE OF MATRIX A ( AND 

R) MCPY 

170 

c 

C - GENERAL 


MCPY 

160 

c 

1 - SYMMETRIC 


MCPY 

190 

c 

2 • DIAGONAL 


MCPY 

200 

c 



MCPY 

210 

c 

REMARKS 


MCPY 

220 

c 

NONE 


MCPY 

230 

c 



MCPY 

240 

c 

SUBROUTINES AND FUNCTION SUBPROGRAMS 

REQUIRED 

MCPY 

250 

c 

LOC 


MCPY 

260 

G 



MCPY 

270 

c 

METHOD 


MCPY 

260 

c 

EACH ELEMENT OF MATRIX A IS MOVED 

TO THE CORRESPONDING 

MCPY 

290 

c 

ELEMENT CF MATRIX R 

. 

MCPY 

300 

c 


1 . , 

MCPY 

310 

c 



> • • •MCPY 

320 

c 



MCPY 

330 


SUBROUTINE MCPY(A,R»N>MrMS) 


MCPY 

340 


DIMENSION A(l)^RU) 


MCPY 

350 

c 



MCPY 

360 

c 

COMPUTE VECTOR LENGTH, IT 


MCPY 

370 

c 



MCPY 

380 


CALL LCC(N,M>1T,N,M9MS) 


MCPY 

390 

c 



MCPY 

400 

c 

COPY MATRIX 


MCPY 

4X0 

c 



MCPY 

420 


DC 1 I*ltIT 


MCPY 

430 


1 Rn)«Am 


MCPY 

440 


RETURN 


MCPY 

450 


END 


MCPY. 

460 . ; 

1 

i 
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ooooonoooonooor»ooooonooor>ooooooooooor»n o o onoooooo 


SUBROUTINE FSTR 
PURPOSE 

CHANGE STORAGE MODE OF A MATRIX 
USAGE 

CALL MSTR(A,R,N,MSA,MSR) 


DESCRIPTION 
A - 
R - 
N - 
MSA 


OF PARAMETERS 
NAME CF INPUT MATRIX 
NAME CF OUTPUT MATRIX 
NUMBER OF ROWS AND COLUMNS IN A AND 
- ONE DIGIT NUMBER FOR STORAGE MODE 

0 " GENERAL 

1 - SYMMETRIC 
7 - DIAGONAL 

SAME AS MSA EXCEPT FOR MATRIX R 


R 

OF 


MATRIX A 


MSR - 


REMARKS 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX A 
MATRIX A MLST BE A SQUARE MATRIX 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
LOC 


METHOD 

MATRIX A 
MSA MSR 


10 

20 

30 

40 

50 

60 

70 

80 

90 


IS RESTRUCTURED TC FORM MATRIX R* 

MATRIX A IS MOVED TO MATRIX R 

THE UPPER TRIANGLE ELEMENTS OF A GENERAL MATRIX 
ARE USED TO FORM A SYMMETRIC MATRIX 
THE DIAGONAL ELEMENTS OF A GENERAL MATRIX ARE USED 
TC FORM A DIAGONAL MATRIX 

A SYMMETRIC MATRIX IS EXPANDED TO FORM A GENERAL 
MATRIX 

MATRIX A IS MOVED TO MATRIX R 

THE DIAGONAL ELEMENTS OF A SYMMETRIC MATRIX ARE 
USED TO FORM A DIAGONAL MATRIX 

A DIAGONAL MATRIX IS EXPANDED BY INSERTING MISSING 
ZERO ELEMENTS TO FORM A GENERAL MATRIX 
A DIAGONAL MATRIX IS EXPANDED BY INSERTING MISSING 
ZERO ELEMENTS TO FORM A SYMMETRIC MATRIX 
MATRIX A IS MOVED TO MATRIX R 


MSTR 
,MSTR 
MSTR 
MSTR 
MSTR 
MSTR 
MSTR 
MSTR 
MSTR 
MSTR IOC 
MSTR 110 
MSTR 120 
MSTR 130 
MSTR 140 
MSTR ISO 
MSTR 160 
MSTR 170 
MSTR 180 
MSTR 190 
MSTR 200 
MSTR 210 
MSTR 220 
MSTR 230 
MSTR 240 
MStR 250 
MSTR 260 
MSTR 270 
MSTR 280 
MSTR 290 
MSTR 300 
MSTR 310 
MSTR 320 
MSTR 330 
MSTR 340 
MSTR 350 
MSTR 360 
MSTR 370: 
MSTR 360 
MSTR 3<90 
MSTR 400 
MSTR 410 
MSTR 420 
MSTR 430 
MSTR 440 
MSTR 450 
MSTR 460 
MSTR 470 
.MSTR 480 
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C ‘ ■ 

SUBROUTINE KSTR ( A,R,N,MSAi MSR) 

DIMENSION AM),R(1} 

C 

DC 20 1«1,N 
DO 20 J*l,N 

IF R IS GENERALf FORM ELEI'ENT 
IF(MSR) 5,10»5 

IF IN L0V«ER TRIANGLE OF SYMMETRIC OR DIAGONAL Rt BYPASS 

5 IF(I-J) 10,10,20 
10 CALL LOC(I,J,IR,N,N,MSR) 

IF IN UPPER AND OFF DIAGONAL OF DIAGONAL R, BYPASS 

IF( IR) 20,20,15 

OTHERWISE, FORM RlltJ) 

15 R( IR)«0.0 

CALL LOC(i,J,lA,N,N, MSA) 

IF THERE IS NC A(I,J), LEAVE R(lfJ) AT 0.0 

IF(IA) 20,20,18 

16 R( IR)>A( U) 

20 CONTINUE 

RETURN 
END ■ 


MSTR 490 
MSTR 900 ! 
MSTR 910 
MSTR 920 
MSTR 930 
MSTR 940 
MSTR 990 
MSTR 940 
MSTR 9T0 
MSTR 9B0 
MSTR 990 
MSTR 600 
MSTR 610 
MSTR 620 
MSTR 630 
MSTR 640 
MSTR 650 
MSTR 660 
MSTR 670 
MSTR 6B0 
MSTR 690 
MiTR 700 
MSTR 710 
MSTR 720 
MSTR 730 
MSTR 740 
MSTR 790 
MSTR 760 
MSTR 770 
MSTR 760 
MSTR 790? 
MSTR 800 
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SUhROUTIivJE loc 


c 

c 

c 

c 

c 

c 

c 


PURPOSE 

C0!^PUTE 


A VPCTOk SUmSCH^PT hop fiN ELEMENT IN A MATRIX OE 


G 

SPECrFT(“0 STORAGE M(U)F 

C 



C 

usage 


C 

CALL 

LUC ( I * J * IR *N * M . MS ) 

C 



C 

Dff;srR I PT T ON OF PARaMETE'RS 

c 

1 

k(jw NUMBER OF element 

C 

J 

COLUMN NUmhER of EILEmFNT 

c 

TR - 

resultant VeCTQP SUHSCPIPT 

c 

N - 

fvjUMHF^ OF ROi^S IN MATRIX 

c 

M — 

number of COLUMNS IN MATRIX 

c 

mS 

ONE DIGIT number for STORAGE 

C 


0 - GFNEhAL 

c 


1 - symmetric 

c 


P - DIAGONAL 

C 



c 

REMARKS 


C 

NONE 



MODE OF MATRIX 


SURPOUT TNES 
NONE 


ANO FUNCTION SUHRROORAMS REQUIRED 


method 

MS = 0 
MS=1 


MS = ? 


SUHSOIRI IS COMPUTED FOR A MATRIX WITH N«“M ELEMENTS 

IN storage (oenfral matrix) 

SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N*<N+1)/? IN 
storage (UPPER TRIANGLE OF SYMMETRIC MATRIX). IF 
ELEMENT IS IN LOWER TRIANGULAR PORTION* SUBSCRIPT IS 
CORRESPONDING ELEMENT IN UPPER TRIANGLE. 
subscript is computed for a matrix WITH N ELEMENTS 
IN STORAGE (DIAGONAL ELEMENTS OF DIAGONAL MATRIX). 

IF element IS NOT ON DIAGONAL (AND THEREFORE NOT IN 
STORAGE)* IR IS SET TO iTERQ. 


10 

?0 

22 

?4 

TO 

32 

36 


SUhROUTINF loc (I * J» ip »N*M»ms) 

TX = I 
JX^J 

IF(MS-l) 10*20*30 
I^X=N* ( JX- 1 ) + I X 
GO TO 36 

IF(IX-JX) 22*24*2A 
IRX=IX+ (JX*JX-JX)/2 
GO TO 36 

IPX=JX+ ( I X«IX-I X) /2 
GO TO 36 
IRX = Q 

IF ( IX-JX) 36*32*36 

IRX=IX 

IR=IRX 

RETURN 

END ‘ 


LOC 

10 

,LOC 

20 

LUC 

30 

LOC 

40 

LOC 

50 

LOC 

60 

LOC 

70 

LOC 

80 

LOC 

90 

LOC 

100 

LOC 

no 

LOC 

120 

LOC 

130 

LOC 

140 

LOC 

150 

LOC 

160 

LOC 

170 

LOC 

180 

LOC 

190 

LUC 

200 

LOC 

210 

LOC 

220 

LOC 

230 

LOC 

240 

LOC 

250 

LOC 

260 

LOC 

270 

LOC 

200 

LOC 

290 

LOC 

300 

LOC 

310 

LOG 

320 

LOC 

330 

LOC 

340 

LOC 

350 

LOC 

360 

LOC 

370 

LOC 

300 

LOC 

390 

LOC 

400 

LOC 

410 

.LOC 

420 

LOC 

430 

LOC 

440 

LOC 

450 

LOC 

460 

LOC 

470 

LOC 

400 

LOC 

490 

LOC 

500 

LOC 

510 

LOC 

520 

LOC 

530 

LOC 

540 

LOC 

550 

LOC 

560 

LOC 

570 

LOC 

500 

LOC 

590 

LOC 

600 

LOC 

610 
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GMSU 
. .GMSU 

10 

?Q 


1 

GMSU 

30 

SL0RCUTINE Gf'SLB 


GMSU 

40 



GMSU 

50 

PURPOSE 

. 

GMSU 

60 

SUBTRACT CNE GENERAL MATRIX FROM 

ANOTHER TO FORM RESULTANT 

GMSU 

70 

MATRIX 

• 

GMSU 

80 



GMSU 

90 

USAGE 


GMSU 

100 

call GMSLB(A,R,R,NfM) 


GMSU 

110 



GMSU 

120 

description CF PARAMETERS 


GMSU 

130 

A - NAME CF first INPUT MATRIX 


GMSU 

140 

B - NAME Cr SECOND INPUT MATRIX 


GMSU 

150 

R - NAME CF OUTPUT MATRIX 


GMSU 

160 

N - NUMBER GF ROWS IN A,B,R 


GMSU 

170 

M ^ NUMBER OF COLUMNS IN A,B,R 


GMSU 

180 



GMSU 

190 

REMARKS 


GMSU 

200 

ALL MATRICES MUST BE STORED AS GENERAL MATRICES 

GMSU 

210 



GMSU 

220 

SUBROUTINES AND FUNCTION SU6PRCGRAMS 

REOUIRFC 

GMSU 

230 

NONE 


GMSU 

240 



GMSU 

250 

METHOD 


GMSU 

260 

MATRIX e ELEMENTS ARE SUBTRACTED 

FROM CORRESPONDING MATRIX 

AGMSU 

27C 

ELEMENTS 


GMSU 

280 



GMSU 

290 



. «GM^U 

300 



GMSU 

310 

SUBROUTINE GMSL B ( A , B t R »N , M ) 


GMSU 

320 

DIMENSION A< 1) ,B( 1) ,R{U 


GMSU 

330 



GMSU 

340 

CALCULATE NUMBER OF ELEMENTS 


GMSU 

350 



GMSU 

360 

NM=N*M 


GMSU 

370 



GMSU 

360 

SUBtRACT MATRICES 


GMSU 

390 



GMSU 

400 

DC 10 1*1, NH 


GMSU 

410 

Rm*A(n-Bm 


GMSU 

420 

RETURN 


GMSU 

430 

END 


GMSU 

440 
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C 


GKPR 10 

G^'PR 20 

GKPR '30 

SLBRCUTINE G^PRD ^MPR AO 

Gf^PR 50 

PURPOSE GNPR.60 

MULTIPLY TUC GENERAL MATRICES TO FORM A RESULTANT GENERAL GMPR 7C 

MATRIX SO 

GMPR 9C 

USAGE GMPR 100 

CALL GMPRC(A,B,RtN,MtL) Ot^PR 110 


GMPR 120 

DESCRIPTION CF PARAMETERS 

A ' NAME CF FIRST INPUT MATRIX 
B - NAME CF SECOND INPUT MATRIX 
R - NAME CF OUTPUT MATRIX 
N “ NUMBER CF ROWS IN A 

M - NUMBER CF COLUMNS IN A AND RCWS IN B 
U - NUMBER CF COLUMNS IN E 


REMARKS 

ALL MATRICES MUST BE STORED AS GENERAL MATRICES GMPR 220 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX A . GMPR 230 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX B GMPR 240 

NUMBER CF columns OF MATRIX A MUST BE EQUAL TO NUMBER OF ROWGMPR 250 
OF MATRIX B , GMPR 260; 

GMPR 270 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED GMPR 280 

NONE ■ GMPR 290 

GMPR 300 

METHOD GMPR 310 

THE M BY L MATRIX B IS PREMULTIPLIED BY THE N BY M MATRIX A GMPR 320 

AND THE RESULT IS STORED IN THE N BY L MATRIX R* GMPR 330 

GMPR 340 

..••.GMPR 3^0 


GMPR 130 
GMPR 140 
GMPR 150 
GMPR 160. 
GMPR 170 
GMPR leo 
GMPR 19C 
GMPR 2CC 


SUBROUT INE GMPR C M , B , R , N » M , L 
DIMENSION A( 1 ) tB ( U ,R( 11 


IR=0 
IK = -M 

CO 10 K=l,L 
IK*IK+M 
CO 10 J=ltN 
IR=1R+1 
JI = J~N 

R( IR )*0 


GMPR 360 
GMPR 37C 
GMPR 380 
GMPR 390 
GMPR 400 
GMPR 410 
GMPR, 420 
GMPR 430 
GMPR 440 
GMPR 450 
GMPR 460 
GMPR 470 
GMPR 480 
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I 


c 


$U8RCUTINE GKTRA 
PURPOSE 

TRANSPOSE A GENERAL p!ATRIX 
USAGE 

CALL GHTRA(A,R,N,K) 

DESCRIPTION CF PARAMETERS 

A - NAME CF MATRIX TO BE TRANSPOSED 
R - NAME CF RESULTANT MATRIX 
N - NUMBER OF ROWS OF A AND COLUMNS OF R 
M - NUMBER OF COLUMNS OF A AND ROWS OF R 

REMARKS 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX A 
MATRICES A AND. R MUST BE STORED AS GENERAL MATRICES 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

TRANSPOSE N BY M MATRIX A TO FORM M B.V N MATRIX R 


SUBROUTINE GMTR A t A ,R »N» M ) 
DIMENSION A(1),RU) 

IRsQ 

DC 1C lelfN 
IJ*I-N 
DC 1C J=lfM 
IJ*U+N 
IR«IR4l 

10 R(IR)*AUJ) 

RETURN 

END 


GMTR 10 
,GMTR 20 
GMTR 30 
GMTR ^0 
GMTR 50 
GMTR 60 
GMTR 70 
GMTR 80 
GMTR 90 
GMTR 100 
GMTR 110 
GMTR 120 F 
GMTR 130 
GMTR lAO 
GMTR 150 
GMTR 160 
GMTR 170 
GMTR 180 
GMTR 190 
GMTR 200 
GMTR 210 
GMTR 220 
GMTR 230 
GMTR 2AC 
GMTR 250 
GMTR 260 
GMTR 270 
.GMTR 280 
GMTR 290 
GMTR 3C0 
GMTR 310 
GMTR 320 
GMTR 330 
GMTR 340 
GMTR 350 
GMTR 360 
GMTR 370 
GMTR 380 
GMTR 390 
GMTR 400 
GMTR 410 




CC 1C 1*1, M 
Jl= JI+N 
ie=ie+i 

IC R nR)=RnH)+A { jn^BI IB) 
RETURN 
END 


G^PR ^9C 
GMPR 500 
Gf^PR 510 
GHPR 520 
G^PR 530 
GMPR 540 
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SIBRCUTINE KPRD 


PURPCSE 

MULTIPLY 


me MATRICES TO FCRM A RESULTANT MATRIX 


USAGE 

CALL MPRC ( A,Bf R,N,N,MSA,MSB,L) 
DESCRIPTICN CF PARAMETERS 


A - 
B - 
R - 
N • 
M - 
MSA 


MSB 

L - 


name CF FIRST INPUT MATRIX 

NAME CF SECOND INPUT MATRIX 

NAME CF OUTPUT MATRIX 

NUMBER CF ROWS IN A AND R 

NUMBER OF COLUMNS IN A AND ROWS IN B 

- CNE DIGIT NUMBER FOR STORAGE MODE OF MATRIX A 

C - GENFRAL 

1 - SYMMETRIC 

2 - CIAGUNAL 

- SAME AS MSA EXCEPT FCR MATRIX B 
NUMBER OF COLUMNS IN B AND R 


REMARKS 

MATRTX k CANNOT BE IN THE SAME LOCATION AS MATRICES A OR 
NUMBER OF COLUMNS OF MATRIX A MUST BE EQUAL TO NUMBER OF 
OF MATRIX B 


SUBROUTINES 

LCC 


AND FUNCTION SUBPROGRAMS REQUIRED 


METHOD 

THE M BY t MATRIX B IS PREMULT 1 PL I ED 
AND THE RESULT IS STORED IN THE N BY 
ROW INTO CCLUMN PRODUCT. 

THE FOLLCWINC TABLE SHOWS THE STORAGE 
MATRIX FCR ALL COMBINATIONS OF INPUT 
A B 

GENERAL 
SYMMETRIC 
DIAGONAL 
GENERAL 
SYMMETRIC 
DIAGONAL 
GENERAL 
SYMMETRIC 
DIAGONAL 


GENERAL 

GENERAL 

GENERAL 

SYMMETRIC 

SYMMETRIC 

SYMMETRIC 

DIAGONAL 

DIAGONAL 

DIAGONAL 


BY THE N BY M MATRI 
L MATRIX R. THIS IS 

MODE OF THE OUTPUT 

matrices 

R 

GENERAL 

GENERAL 

GENERAL 

GENERAL 

GENERAL 

GENERAL 

GENERAL 

GENERAL 

DIAGONAL 


MPRC 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRC 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
ROWMPRD. 
MPRD 
^ MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRD 
MPRO 
MPRD 
MPRD 
MPRO 
MPRD 
MPRD 
MPRO 
MPRO 
MPRO 


B 


X A 
A 


1C 

20 

30 

AO 

5C 

60 

70 

80 

90 

100 

110 

120 

130 

lAO 

150 

160 

170 

180 

190 

200 

21C 

220 

230 

2A0 

250 

260 

270 

280 

290 

300 

31b 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460' 

470 

480 
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SUBRCUTINE MPRC( AfBtRfNiM.MSAff'SB ,U 
CI^ENSICN A( 1) ,e ( 1 ) ,R( 1 ) 

special case for diagonal by diagonal 


f^S*KSA*10+MSB 
IF(^^S-22) 30,lCt30 
10 DC 20 1 = 1, N 
20 R( I)=A( I )«R( I ) 

RETURN 

ALL OTHER CASES 

30 IR=1 

DC 90 K=1,L 
DC 9C J=1,N 
R{ IR )=0 
CG eC 1 = 1,1^ 
lE(KS) AO, 60, AC 
AO CALL LOC( J,I ,IA,N,M,MSA) 
CALL LOC( I ,K, ie,^^,L,H$R) 
IF(IA) 5C,8C,5C 
50 IF ( IR) 7C,«0,7C 
60 IA=N«(I-1)+J 
IB=K*(K-1)+I 

70 RUR)=R( IR) + A(T AI’i'Bnb) 
80 CONTINUE 
90 IR=IR+1 
RE TURN 
END 


MPRC A90 
.NPRD 50C 
HPRC 510 
HPRC 520 
KPRD 530 
KPRD 5A0 
HPRO 550 
NPRC 560 
MPRD 570 
MPRD 580 
MPRD 590 
MPRD 600 
MPRD 610 
MPRO 620 
MPRD 630 
MPRD 6AC 
MPRD 650 
MPKO 660 
MPRD 670 
MPRD 680 
MPRD 690 
MPRD 700 
MPRD 710 
MPRD 720 
MPRO 730 
MPRD 7A0 
MPRD 750 
MPRD 760 
MPRO 770’ 
MPRD 780 
MPRD 790 
MPRO 800 
MPRD 810 
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r 

G 

G 

C 

C 

c 

G 

C 

C 

C 

c 

c 

G 

C 

c 

c 

c 

c 

c 

c 

G 

C 

C 

c 

c 

G 

C 

C 

G 

C 

C 

C 

c 

G 

C 

C 

C 

G 

C 

C 

C 

c 


Su►^r^')UT7^'F CCur 


SUHQOUTP'jF CCUT 


FUiVPdS*’ 

•^EnULTanT 

USAf4F. 

CALL CCUT 


A »«1AThIX 
^ATkICES 


HETwEtN SPtCIFIfcO COLUMNS TO FOHM TWO 




OFSCPIMriiiN 
A 
L 


S 

N 

.M . 

MS 


OF PAwAMf TEPS 
namf of INmUT -^atpix 
COLUMN <)F A TO The LEFT 
PLACE 

name of waTpIX TO HE 

NAME OF MAT^rx TO HE 

numpep of kows tn a 

NUMNEP OF COLUMNS IN A 
- ONE OIHIT NO^iRfcH FOR 

0 - oenfhal 

1 - SYMMF Th T c 

? - i-UAOONAL 


OF which PARTITIONING TAKES 


formed from left PORTION OF A 
formed from right portion of a 


STORAGE MODE OF MATRIX A 


REWARDS 

matrix r cannot me in same 

matrix S cannot he IN SAME 

matrix r Cannot me in same 

matrix R ANI) matrix S AkE 


LOCATION 

LOCATION 

location 


AS 

AS 

AS 


MATRIX 

MATRIX 

matrix 


SUBROUTINES 

LOC 


ALWAYS GENERAL MATRICES 
AND FUNCTION SUhPROGHAMS REQUIRED 


MFTHOf^ 

ELEmEmTS of 
form matrix 

MATRIX A IN 

matrix s of 


matrix a to the left of column L ARE MOVED TO 
R OF N RO'^S aNO L-i COLUMNS. ELEMENTS OF 
COLUMN L AND TO THE RIGHT OF L ARE MOVED TO F( 
N ROWS AND M-L+1 COLUMNS. 


SUBROUTINE CCUT ( A *L ♦ R* S ♦ Nf M* MS ) 
niMENSTON AM)*R(l)»S(n 

TPsO 


CCUT 

10 

. . .CCUT 

20 

CCUT 

30 

CCUT 

40 

CCUT 

50 

XCUT 

60 

CCUT 

70 

CCUT 

80 

CCUT 

9.0 

CCUT 

100 

CCUT 

110 

CCUT 

120 

CCUT 

130 

CCUT 

140 

CCUT 

150 

CCUT 

160 

CCUT 

170 

CCUT 

180 

CCUT 

190 

CCUT 

200 

CCUT 

210 

CCUT 

220 

CCUT 

230 

CCUT 

240 

CCUT 

250 

CCUT 

260 

CCUT 

270 

CCUT 

280 

CCUT 

290 

CCUT 

300 

CCUT 

310 

CCUT 

320 

CCUT 

330 

CCUT 

340 

CCUT 

350 

CCUT 

360 

CCUT 

370 

RMCCUT 

380 

CCUT 

390 

CCUT 

400 

..CCUT 

410 

CCUT 

420 

CCUT 

430 

CCUT 

440 

CCUT 

450 

CCUT 

460 
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1^=0 
nr> 7o 

I >o 7 0 j = 1 , iM 

FIND I. OC^IIOM I ‘^i <’UTHi,JT (WfilKlX ANO SET TO ZERO 

( J-L ) » 1 0 . I 0 

10 TS=IS*1 
S ( T S ) = 0 . 0 
GO TO LID 
?0 17^+1 

( T ^ ) = 0 , 0 

L()CAT7' FLI:^^EN( FOh> any maTkJX STOHAOE N»00F 
30 CALL LOC { T 9 J? IJ*09W«Mb) 

TEST FO^< ZFhO FLF'^LNT In OIAfiONAL MAT»^IX 
T ( IJ ) 4 0*70940 

u>ftfmn*in£ whetheh wight o«' left of l 

40 TF^J-L) '^OtbOtSO 
s f» S { T S ) ^ A ( IJ ) 

GO TO 70 
hO P(Tf<)=AU.j) 

70 CONTINUF 
pFTU^^n 
F N n 


CCUT 470 
CCUT 4ft0 
CCUT 490 
CCUT 500 
CCUT 510 
CCUT 520 
CCUT 530 
CCUT 540 
CCUT 550 
CCUT 560 
CCUT 570 
CCUT 580 
CCUT 590 
CCUT 600 
CCUT 610 
CCUT 620 
CCUT 630 
CCUT 640 
CCUT 650 
CCUT 660 
CCUT 670 
CCUT 680 
CCUT 690 
CCUT 700 
CCUT 710 
CCUT 720 
CCUT 730 
CCUT 740 
CCUT 750 
CCUT 760 
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StBROUTINE >*INV 
PURPOSE 

INVERT A KATRIX 
USAGE 

CALL KINV(AtN«0«L»P) 

DESCRIPTION OF PARAMETERS 

A ^ INPUT MATRIX^DESTROYEC IN COMPUTATION AND REPLACEC EY 
RESULTANT INVERSE* 

N - ORDER OF MATRIX A 
0 - RESULTANT DETERMINANT 
L - WORM VECTOR OF LENGTH N 
M - WORK VECTOR OF LENGTH N 

REMARKS 

MATRIX A MUST BE A GENERAL MATRIX 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

the STANDARD GAUSS- JORDAN METHOD IS USED* THE DETERMINANT 
IS ALSO CALCULATED. A DETERMINANT OF 2ERO INDICATES THAT 
THE MATRIX IS SINGULAR. 


SUBRCUTINE M INV ( A •NtOtL t M ) 

dimension AU)iLm*Mm 


IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESlREOt THE 
C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION 
STATEMENT WHICH FOLLOWS* 

DOUBLE PRECISION A»D»BI6A, HOLD 

THE C MUST ALSO BE REMOVED FRCM DOUBLE PRECISION STATEMENTS 
APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS 
ROUTINE. 

THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO 


MINV 

»MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV. 

MINV 

MINV 

MINV 

MINV 

MINV 

MITIV 

MINV 

MINV 

.MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 

MINV 


10 

20 ' 

30 

40 

50 

60 

70 

80 

90 

100 

no 

120 

130 

140 

130 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

4S0 

4^0 

470 

480 
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CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. 
10 MUST BE CHANGED TO DABS. 


SEARCH FOR LARGEST ELEMENT 


C-1.0 

NK*-N 

DC 80 K»ltN 

NK*NK+N 

UK 

M(K)«K 

KK«NK4^K 

BIGA«A(KK) 

CG 20 J>K«N 
IZ»N*U-1) 

00 20 I*K,N 
IJ*IZ+I 

10 IF( ABS(BIGA)' AeSIAdJU) 15«20»20 
15 BICA»A( IJ ) 

L(K)«I 
M(K»«J 
20 CONTINUE 

INTERCHANGE RCNS 

J*LIK) 

IFU-K) 35f35,25 
25 KI*K-N 

CC 30 I«ltN 

KI*KI+N 

H0LD=^ACKI) 

J!«KI-K+J 
A(KI)»A(jn 
30 AUn -HOLD 

INTERCHANGE CCLUMNS 

35 I»M(K) 

IF(I-K) 45,45,38 
38 JP*N*(I-l) 

CO 40 J«I,N 
JK-NK^J 

HOLD— AUK) 

A(JK)-A(JI) 

40 AUl) -HOLD 


ABS IN STATEMENT MINV 490 

MINV 500 
MINV 510 

..•••MINV 520 

MINV 530 
MINV 540 
MINV 550 
MINV 560 
MINV 570 
MINV 580 
MINV 590 
MINV 600 
MINV 610 
MINV 620 
MINV 630 
MINV 640 
MINV 650 
MINV 660 
MINV 670 
MINV 680 
MINV 690 
MINV 700 
MINV 710 
MINV 720 
MINV 730 
MINV 740 
MINV 750 
MINV 760 
MINV 770 
MINV 780 
MINV 790 
MINV 800 
MINV 810 
MINV 820 
MINV 830 
MINV 840 
MINV 850 
MINV 860 
MINV 870 
MINV 880 
MINV 890 
MINV 900 
MINV 910 
MINV 920 
MINV 930 
MINV 940 
MINV 950 
MINV 960 
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DIVIDE COLUMN BV MINUS PIVOT (VALUE OF PIVOT ELEMENT IS 
CONTAINED IK BIGA) 

A5 IF(BIGA) A8,A6,48 
A6 C»C.O 
RETURN 

40 DC 55 I=1*N 

IF(I-K) 5C»55f5C 
50 IK*NK+I 

A( IK)«A( IK)/(-6IGA) 

55 continue 

REDUCE MATRIX 

DC 65 I«l,N 

IK*NK+I 

H0LD»A(IK) 

IJ*I-N 
DC 65 J*lfN 
1J*IJ4N 

IF(I-K) 6C,65»6C - 

60 IFU-K) 62,65,62 
62 KJ-IJ-I+K 

A(IJ)»HQLC4A(KJ)-^Anj) 

65 continue 

DIVIDE ROW BY PIVOT 

KJ«K-N 
DC 75 J»1,N 
KJ*KJ+N 

IF(J-KI 70,75,70 
70 A(KJ)«A(KJ)/8ICA 
75 CONTINUE 

PRODUCT OF PIVOTS 

0«D»BI6A 

REPLACE PIVCT BY RECIPROCAL 

A(KK)«1.0/eiGA 
80 CONTINUE 

FINAL ROW AND COLUMN INTERCHANGE 

K*N 


MINV 970 
MINV 980 
MINV 990 
MINVIOOO 
MIN VI 010 
MINV1020 
MINV1030 
MINV1040 
MINV1050 
MINV1060 
MINV1070 
MINV1060 
M1NV1090 
MINVllOO 
MINVIXIO 
MINV1120 
MINV1130 
MINV1140 
M1NV1150 
M1NVU60 
MINVU70 
M1NV1180 
MINVU90 
MINV1200 
M1NV1210 
MINV 1220 
MINV1230 
M1NV1240 
M1NV1250 
MINV 1260 
MINV1270 
MINV1280 
MINV1290 
M1NV1300 
MINVUIO 
M1NV1320 
MINV1330 
MINV1340 
MINV1350 
M1NV1360 
MINV1370 
MINV1380 
MINV1390 
MINV1400 
M1NV1410 
MINV1420 
MiNV1430 
MINV 1440 
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ICO K*{K-1) 

IF(K) 150il5C,lC5 
1C5 I*UK) 

IF(I-K) 120»l2CilC8 
108 JC»N*(K-n 
JR = N*( I-l ) 

CO no J*l.N 
JK*JC+J 
hClC»A( JK) 

A(JK)*-A(jn 
110 A( JI ) «HOLO 
120 

IF( J-K) lOOf lOCf 125 
125 KI*K-N 

CC 130 1-ltN 

KI*KI+N 

HCLD»A(KI) 

JI«KI-K+J 

A(KI)»-A(JI> 

130 A(JI) »HOLD 
GO TO 100 
150 RETURN 
ENC 


HINV1450 

RINVIA60 

^^INVIATC 

HINVH80 

F1NV1490 

FINV1500 

FINV1510 

FINV1520 

FINV1530 

FINV15A0 

HINV1550 

FINV1560 

FINV1570 

F1NV1580 

F1KV1590 

F1NV160Q 

HINV1610 

FINV1620 

KINV163C 

FINV16A0 

FINV1650 

FINV1660 

FINV1670 

F1NV1680 
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SINV 10 

SINV 20 

SINV 30 

St'BRCUTlNE SINV SINV 40 

SINV 5C 

PURPOSE SINV 60 

INVERT A GIVEN SVPMETKIC POSITIVE DEFINITE MATRIX SINV 70 

SINV 80 

USAGE SINV 90 

CALL SINV (A,N,EPS» lER) SINV 100 

SINV no 

DESCRIPTION CF PARAMETERS SINV 12C 

A - LPPER triangular PART OF THE GIVEN SYMMETRIC SINV 13C 

PCSITIVE DEFINITE N BY N COEFFICIENT MATRIX. SINV 140' 

ON RETURN A CCNTAINS THE RESULTANT UPPER SINV 150 

TRIANGULAR MATRIX. SINV 160 

N - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX. SINV 170 

EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE SINV 180 

TCLERANCE FOR TEST ON LOSS OF SIGNIFICANCE. - SINV 190 

lER - RESULTING ERROR PARAMETER CODED AS FOLLOWS SINV 200 

IER*0 - NO ERROR SINV 210 

I6R = -1 - NO RESULT BECAUSE OF WRONG INPUT PARAME- SINV 220 .. 

TFR N CR BECAUSE SOME RADICANO IS NON- SINV 230 

POSITIVE (MATRIX A IS NOT POSITIVE SINV 240 

DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI- SINV 250 

. FICANCE ) SINV 260 

IER»K - WARNING WHICH INDICATES LOSS OF SIGNIFI- SINV 270 
GANCE. THE RADICAND FORMED AT FACTORIZA- SINV 280 
TION STEP K+1 WAS STILL POSITIVE BUT Nd SINV 290 
LONGER GREATER THAN ABS ( EPS*A( K+l,K+l )) . SINV 300 

SINV 310. 

REMARKS SINV 320 

THE UPPER TRIANGULAR PART CF GIVEN MATRIX IS ASSUMED TO BE SINV 330 
STORED COLUMNWISE IN N*(N+l)/2 SUCCESSIVE STORAGE LOCAT IONS . S INV 340 
IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU- SINV 3$0 
LAR MATRIX IS STORED COLUMNWISE TOO. SINV 360- 

THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL SINV 370 
' calculated RADICANOS are PCSITIVE. SINV 380 

SINV 39C 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED SINV 400 ■ 

MFSC SINV 410 

SINV 420 

METHOD SINV 430 

SOLUTION IS DONE USING THE FACTORIZATION BY SUBROUTINE MFSO.SINV 440 • 

SINV A50 

.v...... .SINV 460 

SINV 470 : 

SUBROUTINE S INV ( A ,N , EPS ♦ I ER ) SINV 480 
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c 

C 

DIMENSION A(l) 

DQL 16 LE PRECISICN DIN, WORK 

FACTORIZc GIVEN MATRIX BY MEANS OF SUBROUTINE MFSD 
A = TRANSPCSem * T 
CALL MFSOI A,N,EPS , lER ) 

IF( lER ) 9,1» 1 

INVERT UPPER TRIANGULAR MATRIX T 
PREPARE INVBPSION-LCCP 

1 IPIV = N’«'(N^-1 ) /2 
IND=IPIV 

INITIALIZE INVERSION-LOOP 
CO 6 1 = 1, N 

CIN=1*D0/CBLE( A ( IPI V) I 
A( IPIV) =CIN 
MIN = N 
KENC=I-1 
LANF=N-KEND 
IFIKENO) 5,5,2 

2 J=INn 
C 

C INITIALIZE RGVv-LOOP 

DO A K=l,KENO 
WCRK=0,CO 
MIN=MIN-1 
LHCR=IPIV 


SINV A90 
SINV 500 
SINV 510 
SINV 520 
SINV 530 
SINV 5A0 
SINV 550 
SINV 560 
SINV 570 
SINV 580 
SINV 590 
SINV 600 
SINV 610 
SINV 620. 
SINV 630 
SINV 6A0 
SINV 650 
SINV 660 
SINV 670 
SINV 680 
SINV 690 
SINV 700 
SINV 710 
SINV 720 
SINV 730 
SINV 7A0 
SINV 750 
SINV 760 
SINV 770 
SINV 780 
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SUS^^OUTINt '-^FSn 


MFSO 

MFSO 


PUWt>OSF 

FACTOP 


fttVFN SYM’^fcTPIC POSITIVE DEFINITE MATRIX 


usapf 

CALL MFSO ( a •M*EPSt lER) 

ofscptpttom Of pawameteps 

t - UPPfR triangular part of the given symmetric 
POSITIVE definite N BY N COEFFICIENt MATRIX. 
ON *^ETURN A CONTAINS Th£ RESULTANT UPPER 

triangular matrix. 

N - The number of rows (COLUMNS) IN GIVEN MATRIX, 

frps - AN INPUT constant WHICH IS USED AS RELATIVE 

TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE. 
TEP - RFSULTINCi F>ROR PARAMETER COOED AS FOLLOWS 


IER=fl - NO ERROR MFSD 

IFR=-1 - NO result because of wrong INPUT PARAME- MFSO 
TER N OR BECAUSE SOME RADICANO IS NON- MFSD 
POSITIVE (MATRIX A IS NOT POSITIVE MFSD 

DEFINITE* POSSIBLY DUE TO LOSS OF SIGNI- MFSD 
FICAnCE) MFSO 

IFP=K - WARNING WHICH INDICATES LOSS OF SIGNIFI- MFSD 
CANCE. THE RAOICAND FORMED AT FACTORIZA- MFSO 
TION STEP K*-l WAS STILL POSITIVE BUT NO MFSO 
LONGER GREATER THAN A8S (EPS*A ( 1 * K> 1 ) ) . MFSO 

' MFSO 

REMARKS MFSO 

The UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE MFSO 
STORED COLUMNWISE IN N*(N4l)/2 SUCCESSIVE STORAGE LOCATIONS.MFSO 
IN the same storage locations The RESULTING UPPER TRIANGU- MFSO 
LAP matrix is stored columnwise too. MFSO 

THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL MFSD 
CALCULATED RADICAND5 ARE POSITIVE. MFSO 

the PRODUCT OF RETURNED DIAGONAL TERMS IS E(.)UAL TO THE MFSO 

S(5UARF-ROOT OF THE DETERMINANT OF THE GIVEN MATRIX. MFSO 

MFSO 

SUBROUTINE'S AND FUNCTION SUBPROGRAMS REQUIRED MFSO 

NONE MFSO 

MFSD 

method MFSO 

SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY. MFSD 

The given matrix is represented as product of two TRIANGULARMFSD 

MATRICES* WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF - MFSO 

The returned right hand factor. MFSO 

MFSD 

•••••MFSD 


IFP=K - 


SUBROUTINE mFSD ( A * N*EPS * lER ) 


dimension A(I) 

DOiiBLE PRECISION OPIV*DSUM 

TESrON WRONG INPUT PARAMETER N 
IP(N-l) 
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1 


C lNITlALI?t* DI AAONAL-LOOP 

KOIV=0 
DO 11 K=1.M 

KPTV=KPIV+K 
TMO=KPJV 
Lfr:''JO = K-l 

c 

t CALCULATF TOLFPANCt 

TOL=AHS(FPS*A(KPrv) ) 

stapt factoki/atiom-loop oven k 
no 11 T=K*'I 
osup=o,no 
TF(LENO) ?•4♦^ 

STAPT T^i^^FP LOOP 
P DO 3 L=ULFNO 
LANF = Kt»I V-L 
Ul»^0=INO-L 

3 nSUM=DSUM+DPLF ( A { L ANF ) *A ( L J WD ) > 

END OF INNfcP LOOP 

TNANSFOPm tLF*^ENT A(INO) 

4 DSUM = OFiLF (A ( 1A»0) )-OSUM 
IF(I-K) 10*5»10 

TEST FOP negative. PIVOT ELE*-£NT 

5 IF (SNGL (OSUM) -TOL ) 6*6*9 
h IF{OSUP) 12*1? *7 

7 IF(IER) 8*8*9 

8 IFP=K-1 

COMPUTE PIVOT ELEMENT 

9 DPT V=DSORT (DSUM) 

A(KPIV)=DPIV 
DPIV=1.00/OPIV 
GO TO 11 

CALCULATE TERMS IN POw 

10 A (IND) =OSUM»OPIV 

11 IND=IND+I 

END OF DIAGONAL-LOOP 
RETURN 

12 IFP=-1 
RETURN 
END 


MFSD 610 
MFSD 620 
MFSO 630 
MFSD 640 
MFSD 650 
MFSD 660 
MFSO 670 
MFSD 680 
MFSD 690 
MFSD 700 
MFSD 710 
MFSD 720 

Th row MFSD 730 

MFSD 740 
MFSD 750 
MFSD 760 
MFSD 770 
MFSD 780 
MFSD 790 
MFSO 800 
MFSD 810 
MFSO 820 
MFSD 830 
MFSD 840 
MFSD 850 
MFSD 860 
MFSO 870 
MFSD 880 

AND FOR LOSS OF SIGNIFICANCE MFSD 890 

MFSD 900 
MFSD 910 
’ MFSD 920 
MFSD 930 
MFSO 940 
MFSD 950 
MFSO 960 
MFSD 970 
MFSO 980 
MFSO 990 
MFSOIOOO 
MFSOlOlO 
MFSO1020 
MFSD1030 
MFSD1040 
MFSO1050 
MFSD1060 
MFSC)1070 
MFSO1080 
MFSO1090 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUPPOUTINE mXO(JT 
PUPPf)SE 

PPOOUCE'^ OUTPUT 
LOGICAL UNIT b 


LISTING OF ANY SIZED ARRAY ON 


USAGE 

CALL MXOUT( ICOOE*A*N*M*MS*LINS*IPOS*ISP) 

OESCPIPTIOn of PAPAMETep'; 

rcOOE- INPUT COUE NUMBER TO RE PRINTED ON EACH OUTPUT PAGE 
A-NAMF OF OUTPUT MATRIX 

n-numhep of rows in a 

M-NUMRER OF COLUMNS IN A 
MS-STORAGE MODE OF A WHEPE mS= 

0- GF^^ERAL 

1- SYmmeTkIC 
^-0 1 AGONAL 

L ins-number OF PRINT LINES ON THE PAGE (USUALLY 60) 
IPOS-NUmBER of print positions across the page (USUALLY 132) 
ISP-LINE SPACING COf)E« 1 FOR SINGLE SPACE* 2 FOR DOUBLE 
SPACF 

REMARKS 

THIS SUBROUTINE HAS BEEN MOOIFIED BY M*MUTTON ON 11/27/71 
TO REDUCE THE AMOUNT OF EXTRA PRINTOUT* TO RETURN THE 
SUBROUTINE TO ITS ORIGINAL FORM MODIFY CAROS 280*480*590* 
660*900*920 ACCORDING TO THE SSP MANUAL AND REMOVE CAROS ' 
281-2.84*691. 


SUBROUTINES 

LOC 


AND FUNCTION SUBPROGRAMS REOUIHEO 


.METHOD 

This subroutine creates a standard output listing of any 

SIZED ARRAY WITH ANY STORAGE MODE. EACH PA(iE IS HEADED WITH 
■ The CODE NUMrtEP*OIMENSIONS AND STORAGE MODE OF THE ARRAY. 
EACH COLUMN AND WOW IS ALSO HEADED WITH ITS RESPECTIVE 
mummer. 


subroutine MXOUT (IC0DE*A*N*M*MS*LINS*IP0S*ISP) 
dimension a (1 ) *B(8) 

1 FORMAT(iHl *.BX* 7HMATRIX *I6*6X*I3*SH R0WS.*6X* 13* 8H COLUMNS* 
18X*13HST0RAGE MODE * 11 * BX * 6HPA(;E *12*/) 

? FORMAT ( 12X*8HCOLUmm * 7 ( 3X U 3 * 1 OX ) ) 

3 format flH ) 

4 FORMATdH *4X*7(F16.6) ) 

5 FORMAT (1H0*4X*7(EI6. 6) ) 

4 FORMATdH *14X*6(E1B.B) ) 

5 format (IHO* 14X*6(F18.«)) 


Jsl 


WRITE HEADING 
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MXOU 

10 

• MXOU 

20 

MXOU 

30 

MXOU 

40 

MXOU 

50 

MXOU 

60 

MXOU 

70 

MXOU 

80 

MXOU 

90 

MXOU 

100 

MXOU 

110 

MXOU 

120 

MXOU 

130 

MXOU 

140 

MXOU 

150 

MXOU 

160 

MXOU 

170 

MXOU 

180 

MXOU 

190 

MXOU 

200 

MXOU 

210 

MXOU 

220 

1 MXOU 

230 

MXOU 

240 

MXOU 

250 

MXOU 

260 

MXOU 

270 

MXOU 

280 

MXOU 

281 

MXOU 

282 

MXOU 

283 

MXOU 

284 

MXOU 

290 

MXOU 

300 

MXOU 

310 

MXOU 

320 

MXOU 

330 

MXOU 

340 

MXOU 

350 

MXOU 

360 

MXOU 

370 

MXOU 

380 

MXOU 

390 

MXOU 

400 

MXOU 

410 

MXOU 

420 

MXOU. 430 

MXOU 

440 

MXOU 

450 

MXOU 

460 

MXOU 

470 

MXOU 

480 

MXOU 

490 

MXOU 

480 

MXOU 

490 

MXOU 

500 

MXOU 

510 

MXOU 

520 

MXOU 

530 

MXOU 

540 
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jPOS/1^-1 
LfMO= (I. INS/ TSP) 

IPaGt=l 
10 lstwt=i 

C ?0 (Gt 1 ) TCOUE »N*m,ms* IPaGE 

?o continue 

JNT= j+NFNn-l 
IPAGt =TPAGE+i 
IF f JMT- m) 33*33«3? 

32 JNT=M 

33 CONTINUE 

C WPtTE (G» 2) ( JCUP* JCUP=Jf JNT) 

IF(ISP-l) 3S,3S*40 
3S VIPTTF(G«3) 

40 LrFNO=LSTPT+LFNO-l 
no 80 L=LSTPT»LTEMn 

FOPM OUTPUT WOW LINF 

DO 5S K = ^,^4FNfi 
KK = K 

JT = J+K-1 

CALL LOC (L ♦ JT ♦ I JNT*N>m«MS) 

P(K) =0.0 

IF (IJnT)50*50*45 
45 P(K)=ft(IJNT) 

50 CONTINUE 

CHFCK IF LAST COLUMN. . IF YES GO TO 60 

IF ( JT-M) S5*6(U60 
55 rowTiNUF 

END OF LINE* NOW WWiTE 

60 IF (ISP-1) 65t65*70 
6S WPITE(8*4) (P ( JW) f JW=1 tKK) 

GO TO 75 

70 WPTTE(6,S> (h ( J to) . J w=1 »KK) 

IE End of poi*(s*go chec^ columns 


.75 IF (N-L ) 85*85*80 
80 CONTINUE 

ENO or PAGE* NOW CHECK FOR MORE OUTPUT 

LSTRT=LsTPT+LtND 
GO TO 20 

END OF COLUMNS* THEN RETURN 

85 IF( JT-M) 90*95*95 
90 J=JT+1 
GO TO 10 
95 return 
END 


MXOU 
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MXOU 
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MXOU 
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730 

MXOU 

740 
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790 
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880 
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MXOU 
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970 
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980 
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990 
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APPENDIX D PROCEDURE FOR ELIMINATING CONSTRAINT 
EQUATIONS IN TRIM PROBLEM 

For linear dynamics and a quadrafic performance criterion the trim problem can be 
written in the form 


0 = a + B6 
r = 1/2 6'R6 


0 ) 

( 2 ) 


with 


a = constant vector of dimension n 
6 = control vector of dimension n 

B = nxm coefficient matrix 
R = mxm positive definite weighting matrix 

The obfective is to find the set of control angles 6 that satisfy (1) and minimize (2), The 
trim solution is given by 

1 

( 3 ) 


6 = - B^a 


# -1 -1 -1 
= R B'(BR B') ’ 

The mxn matrix B is a right inverse of B , i«e., BB = j. 


( 4 ) 


Consider the new trim problem that results from eliminating k of the n equality 
constraints. Suppose that the first k constraint equations in (1) are to be disregarded. The 
problem can always be written in this form by reordering the equations if necessary. Partitioning 
(1) gives that the new trim problem is 


0 = a^ + 6^5^ 


r = 1/2 6' R6 
n n 


(5) 

( 6 ) 


where 


o = 



i 

°1 


°2 


u J 




r- . ’i 


1 

k 

^i 


! k 

B = 



1 <7) 

n-k 

l»2 


■n-k • 

1 


U J 
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The solution to the new trim problem is 


6^ = -R''b^(B2R''b^"'q 


( 8 ) 


The following question is of interest: Without starting the problem over again, 
is it possible to compute 6^ using the solution for 6 ? The answer is affirmative and 
a procedure for computing 6 is developed below. 


From the partitioning (7) of the B matrix 


r''b' 


= _r-'b' i r-’b' ] 


br"’b = 


Taking the inverse of (10) results in 


b,r"’b' I b^r"'b^ 


b^r^’b' [BjR"'^ 


n-k 


( 9 ) 


( 10 ) 


where 




’Qi 

1 

1 

1 

1 

Q 


1 

1 

1 

Q 


Q, = E 


-1 


n-k 


(H) 


and 


= -E''b,R‘’b^(B2R''b^)’’ 

Q3 = (B2R''b^)''b2R’'b'E''b,R''b^(B2R'’b^'’ +(B2R''bP 


E= B,R"'b' -B,R‘’b^(B2R'’bP‘'b2R"'b' 


Premultiplying (11) by (9) yields the right inverse of the B matrix in partitioned form 
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where 


I 

I 


L' 1 “2J 

B* = [x -R‘’b^(B2R‘'b2)"'b2]r''b'e'' (12) 

B* = -[x.r''b^(B2R'’b^)'’b2]r'V,E''b,R‘V2(B2R‘'bP'' + R''b^(B2R*'b^)'’ 


Substituting (12) into (3) and using (7) gives that 

S- B*n -b’^^o 
6 - -B,a, BjOj 


(13) 


If we substitute * 

a, = BjR''B^(B 2 R"’Bp'’a 2 (14) 

into (13) then from (8) and (12) it follows that 

6=6 (15) 

n 

This result states that if the first k elements in the vector a are replaced by the values 
computed from (14) then the solution to the original trim problem becomes the solution to 
the new trim problem created by eliminating the first k constraint equations. 


It is apparent from comparing (11) to (14) that (14) can be replaced by 

0, = -Q, QjOj (16) ' 

this is a more useful equation for computing the new volue of a^ since ond Q 2 
are submatrices of a matrix computed in the solution of. the original problem. 

To summarize, the steps for computing 6 are as follows: 

n 

1) Start the computation of 6 using (3) and (4) in the usual way. 

2) After computing (BR*^B')”\ form the submatrices and according to (11). 

3) Replace subvector a ^ in a by the value computed from (16). 

4) Continue the computation of 6 in the usual way. The result will be 5=6^ . 

The above procedure for computing 6^ does not offer any particular advantage over 
using (8) if the calculations are to be done by hand. If a computer program, on the other hand, 
has been developed to compute 6 then the above procedure minimizes the amount of program 
modification required to compute 5^ , 
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APPENDIX E VERIFtCATIQN OF TRIMS PROGRAM 

Laferal trim of the Space Shuttle is an example of the linear trim problem. The linear 
trim problem is to find the control deflections 6 satisfying the equality constraints 

a + B6 = 0 

and minimizing 

J = 1/2 6'R6 

The solution is 

# 

6 = - B a 

The problem of Space Shuttle trim in roll and yaw (two constraint equations) using the 
following four control deflections: 

• yaw deflection of orbiter engine 1 

e yaw deflection of orbiter engines 2 and 3 

• pitch deflection of orbiter engine 2 

(negative of the pitch deflection of orbiter engine 3) 

• rudder deflection 

was solved at MSFC. The control deflection angles vs flight time for the case when the R 
matrix is 

R = Diag[0.49, 0.49, 0.49, 1.00] 

and the bias torques due to misalignments are 

roll torque = 0.87 x 10^ N*m 

yaw torque = 3.02 x 10^ N*m 

are plotted in Figure El. 

The solution to (supposedly) the same trim problem was also computed using the TRIMS 
program as a check of the program. The resulting plot of control deflection angles vs flight 
time is shown in Figure E2, The TRIMS computation was repeated except without the dorsal 
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fins and the trim solution is plotted in Figure E3, 

The results in Figures E2 and E3 computed by TRIMS do not agree with the results in 
Figure El obtained by MSFC. A comparison of the results does not indicate the reason for 
the difference. The computation of 6 land from a , B , and R in the TRIMS program 
checked against hand calculations. Most likely, the area of difficulty is in the computation 
of the vector a and matrix B from the equotions of motion. 



CONTROL DEFLECTION, 8 (deg) 


Figure El Trim Sol tg ion Computed at MSFC 

□ = YAW DEFLECTION ENGINE 1 V= PITCH DEFLECTION ENGINE 3, AND 

(NEGATIVE OF PITCH DEFLECTION ENGINES) 

0=YAW DEFLECTION ENGINE 2 8 3 0= RUOOER DEFLECTION 



control deflection, 6 (deg) 


t 


■ ' 


Figure E 2 Control Deflections vs Flight Time for Spoce Shuttle Trim in Roll ond Yaw with Addition of Dorsol Fins 



I 

> 


185 


control deflection, 6 (deg) 
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